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Abstract 

Systems  which  vary  significantly  over  an  operating  envelope,  such  as  fighter 
aircraft,  generally  cannot  be  controlled  by  a  single  linear  time-invariant  controller. 
As  a  result,  gain-scheduling  methods  are  employed  to  design  control  laws  which  can 
provide  the  desired  performance.  This  thesis  examines  a  relatively  new  approach 
to  gain-scheduling,  in  which  the  varying  controller  is  designed  from  the  outset  to 
guarantee  robust  performance,  thereby  avoiding  the  disadvantages  of  point  designs. 
Specifically,  the  parameter-varying  (LPV)  aircraft  model  is  linearized  using  linear 
fractional  transformations  (LFT’s),  and  the  resulting  control  problem  is  character¬ 
ized  as  the  solution  to  a  set  of  four  linear  matrix  inequalities  (LMI’s).  The  supporting 
theory  is  reviewed  and  two  pitch-rate  controllers  are  designed;  one  for  the  full  lon¬ 
gitudinal  aircraft  model,  and  another  for  the  short  period  model.  It  is  found  that, 
even  though  the  varying  controllers  are  quite  conservative,  they  can  guarantee  better 
robust  performance  over  a  large  portion  of  an  operating  envelope  when  compared  to 
time-invariant  //-synthesis  controllers. 
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1.  Introduction 

In  the  past  decade,  much  of  modern  control  theory  has  focused  on  determining 
the  ideal  controller  to  optimize  the  performance  of  some  given  system.  This  usually 
requires  simplifying  a  complex  real  world  system  into  a  simpler  workable  mathemat¬ 
ical  model.  This  introduces  many  uncertainties  due  to,  for  example,  neglected  plant 
dynamics,  nonlinearities,  or  inaccurate  knowledge  of  the  model  parameters  or  how 
they  may  vary  during  operation.  To  compensate,  the  designer  must  then  account 
for  these  uncertainties  in  his  system  model.  The  challenge  becomes  one  of  designing 
a  controller  which  can  obtain  as  much  performance  as  possible  in  the  presence  of 
these  modelled  uncertainties.  Inevitably,  there  is  a  tradeoff  between  performance 
and  robustness  and  the  designer  must  then  use  his  skill  to  design  a  controller  which 
can  satisfy  both. 

Generally,  if  the  system  does  not  vary  too  significantly,  it  can  be  modelled  fairly 
accurately  as  a  linear  time-invariant  (LTI)  plant  with  some  small  uncertainties;  in 
such  cases,  standard  H^o  and  //-synthesis  theory  can  be  used  to  find  a  controller 
which  satisfies  both  the  performance  and  robustness  objectives.  Unfortunately,  if 
the  system  varies  considerably  during  operation,  a  satisfactory  level  of  performance 
often  cannot  be  achieved  without  sacrificing  robustness,  and  vice  versa.  Under  such 
demanding  conditions,  designers  usually  adopt  a  gain-scheduling  approach. 

Simply  put,  gain-scheduling  is  a  method  of  changing  the  controller  depending 
upon  the  operating  conditions.  There  are  two  general  schools  of  thought  on  how  to 
generate  these  time- varying  controllers.  The  first  and  more  conventional  approach 
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involves  determining  robustly  performing  controllers  for  several  point  designs,  and 
then  designing  a  method  of  scheduling  the  individual  controllers  such  that  robust 
performance  is  maintained  at  each  design  point.  The  resulting  schedule  can  be  as 
simple  as  a  look-up  table,  or  can  involve  designing  a  mathematical  control  law  relat¬ 
ing  or  approximating  the  individual  controllers  (see,  for  example,  [Bla91],  [AW89], 
[LR93],  and  [SC92]).  While  these  types  of  methods  have  the  advantage  of  designing 
optimal  controllers  for  actual  points  in  the  operating  envelope  using  conventional 
design  methods,  they  also  have  significant  disadvantages.  Foremost,  the  designer 
needs  to  design  several  robust  controllers  to  adequately  cover  the  operating  enve¬ 
lope;  in  doing  so,  he  fixes  the  plant  to  specific  operating  points  and  neglects  its 
time-varying  nature.  As  a  result,  robust  performance  can  no  longer  be  guaranteed 
throughout  the  operating  envelope.  In  addition,  depending  on  the  order  of  the  con¬ 
trollers,  the  resulting  scheduling  problem  is  often  quite  complex.  The  individual 
controllers  can  be  simplified,  but  this  often  comes  at  the  expense  of  performance 
unless  the  parameter  variations  are  moderate  to  begin  with.  In  addition,  robust 
stability  throughout  the  operating  envelope  can  only  be  guaranteed  under  certain 
specific  conditions  [SA91a],[SA91b].  In  fact,  much  of  the  gain-scheduling  research 
done  in  the  past  has  focussed  on  addressing  these  issues. 

A  new  school  of  thought  has  emerged  in  recent  years  that  essentially  reverses 
the  classical  gain-scheduling  problem.  Instead  of  designing  several  point  controllers 
first  and  then  trying  to  determine  a  relationship  between  the  controllers  and  the 
changing  operating  parameters,  these  methods  focus  on  designing  the  relationship 
between  the  gain-scheduled  controller  and  the  changing  parameters  from  the  outset. 
Once  this  relationship  is  known,  appropriate  controllers  for  all  operating  conditions 
can  then  be  determined  from  values  of  the  parameters.  In  short,  these  approaches 
directly  design  time-varying  controllers.  An  obvious  advantage  is  that  the  time- 
consuming  process  of  designing  point  controllers  and  then  a  schedule  is  completely 
avoided;  in  addition,  robust  performance  can  be  guaranteed  over  the  entire  design 
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envelope.  However,  these  methods  all  depend  on  determining  relationships  between 
a  linear  time-varying  plant  and  the  varying  parameters.  While  the  relationship 
between  the  plant  and  the  parameters  can  be  expressed  in  several  different  ways, 
to  date  this  invariably  introduces  conservatism  at  some  point  in  the  process  which 
degrades  the  achievable  level  of  robust  performance. 

Specifically,  the  method  developed  by  Packard,  Apkarian,  and  Cabinet  [Pac94], 
[AG95]  (on  which  this  thesis  is  based)  models  linear  parameter- varying  (LPV)  plants 
and  controllers  as,  respectively,  LTI  plants  and  LTI  controllers  linearly  dependent 
upon  a  separate  varying  parameter  matrix.  This  parameter  block  is  then  treated 
as  an  uncertainty  and  the  resulting  Hoo  design  problem  is  solved  for  the  desired 
controller  using  Small  Gain  theory.  Unfortunately,  due  to  the  Small  Gain  Theorem, 
substantial  conservatism  is  introduced  because  the  resulting  controller  is  designed 
to  address  not  only  real  parameter  values,  but  also  complex  parameter  values  which 
will  not  exist  for  most  physical  systems. 

Another  method,  developed  by  Becker  [BP94],  avoids  breaking  up  the  LPV 
plant;  instead,  the  LPV  state-space  of  the  plant  is  kept  within  the  constraints  of  the 
optimization  problem.  This  avoids  the  conservativeness  of  the  Small  Gain  Theorem, 
but  requires  solving  an  infinite  number  of  constraints  since  the  parameters  are  con¬ 
tinuous.  One  way  to  resolve  this  is  to  discretize  the  parameters  and  formulate  the 
constraints  at  each  combination  of  parameter  values.  The  resulting  stack  of  finite 
constraints  can  then  be  solved  for  the  desired  controller.  In  particular,  if  the  state- 
space  matrices  are  afSne  functions  of  the  parameters,  then  it  has  been  shown  that 
only  the  constraints  corresponding  to  the  vertices  of  the  operating  envelope  need  to 
be  solved.  This  entire  approach,  however,  assumes  that  all  convex  combinations  of 
the  LPV  state-space  matrices  exist  in  the  modelled  system,  which  is  not  necessarily 
true  and  again  introduces  substantial  conservatism. 

Finally,  a  third  approach  introduced  by  Wu  [WYPB94]  extends  the  previous 
one  to  account  for  the  rates  of  parameter  variations.  The  first  two  approaches  do 


1-3 


not  place  any  bounds  on  the  rates  of  the  parameter  variations,  and  the  controllers 
are  then  designed  to  address  any  rate  of  change,  no  matter  how  unrealistic.  In 
Wu’s  approach,  bounds  on  the  parameter  rates  of  change  are  included  and  form 
part  of  the  constraints  to  be  solved.  Unfortunately,  this  method  also  results  in  an 
infinite  number  of  constraints;  in  addition,  the  set  of  matrix  functions  from  which 
the  controller  is  extracted  is  also  infinite.  Thus,  while  it  does  address  parameter 
rates  of  change,  it  suffers  from  all  the  disadvantages  of  Becker’s  approach  and  its 
numerical  tractability  is  even  more  difficult. 

1.1  Objectives 

The  main  objective  of  this  thesis  is  to  apply  some  of  the  latest  techniques  in 
gain-scheduling  to  a  realistic  problem.  Specifically,  an  extension  of  the  method  of 
Packard,  Apkarian,  and  Gahinet  for  uncertain  LPV  systems  [AG95],  [SLBB96]  will 
be  used  to  design  a  robustly  performing  pitch  rate  controller  for  an  F-18  Superma- 
neuverable  fighter  aircraft.  This  extension  incorporates  the  advantages  of  ^-synthesis 
to  handle  remaining  uncertainties  and  add  even  more  flexibility  to  the  design.  Using 
these  techniques,  controllers  for  both  the  full  longitudinal  state-space  aircraft  model 
and  the  short  period  aircraft  model  will  be  designed,  simulated,  and  then  compared. 

In  support  of  this  thesis,  most  of  the  required  theory  has  been  included.  A 
secondary  objective  was  to  become  familiar  with  the  underlying  theory  of  linear 
matrix  inequalities  (LMI’s),  which  are  now  often  used  to  express  the  constraints  of 
complex  modern  control  problems.  This  work  represents  a  very  practical  use  of  the 
newly  released  MATLAB  LMI  Control  Toolbox  [MAT],[GNLC92]. 

1.2  Outline 

Following  this  introduction,  the  remainder  of  the  thesis  is  divided  into  five 
chapters.  Chapter  II  introduces  some  of  the  mathematical  and  control  theory  fun¬ 
damentals  essential  to  the  thesis.  While  the  reader  is  expected  to  be  familiar  with 
general  Hoo  theory,  this  chapter  reviews  the  concepts  of  linear  fractional  transfor- 
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mations  (LFT’s),  LMI’s,  and  //-synthesis  which  will  be  used  extensively  in  later 
chapters. 

Chapter  III  describes  the  gain-scheduling  theory  for  LPV  systems  introduced 
by  Packard,  Apkarian,  and  Cabinet  [Pac94],  [AG95].  It  also  describes  the  extension 
for  uncertain  LPV  systems  mentioned  previously  [SLBB96].  All  of  the  constraints 
to  solve  the  resulting  problems  are  formulated  as  LMFs  which  must  then  be  solved 
to  obtain  the  LPV  controller. 

Chapter  IV  details  the  implementation  of  the  theory  to  solve  for  the  gain- 
scheduled  pitch  rate  controller  for  the  F-18  aircraft  under  study.  The  entire  setup  of 
the  problem  is  described,  beginning  with  how  to  obtain  an  LPV  plant  representation 
and  how  to  then  convert  the  LPV  plant  into  an  LTI  plant  linearly  dependent  upon  a 
varying  parameter  block.  The  design  model  is  then  developed;  for  this  particular  de¬ 
sign,  a  model-matching  framework  is  adopted  with  a  2-degree-of-freedom  controller. 
The  selection  of  weights  and  design  envelopes  are  then  discussed,  along  with  some 
practical  implementation  details  which  may  not  be  readily  apparent. 

Chapter  V  presents  the  numerical  results  and  simulations.  The  resulting  con¬ 
trollers  are  described  and  then  simulated  to  examine  their  performance  throughout 
and  beyond  their  design  envelopes.  Various  types  of  simulations  are  performed  and 
evaluated;  as  well,  the  time  responses  of  LTI  //-synthesis  controllers  is  included  for 
comparison  purposes. 

Finally,  Chapter  VI  summarizes  the  results,  discusses  some  of  the  problems, 
and  recommends  some  areas  for  further  study. 

Several  appendices  provide  further  information.  Appendix  A  lists  the  aircraft 
flight  conditions  used  as  data  in  order  to  establish  a  relationship  between  the  aircraft 
model  and  the  varying  parameters;  Appendix  B  illustrates  the  resulting  relationships. 
As  well,  most  of  the  MATLAB  programs  for  the  implementation  of  the  control 


1-5 


problem  have  been  provided  in  Appendix  C.  Finally,  the  simulation  programs  are 
included  in  Appendix  D. 
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II.  Preliminaries:  LFT’s,  LMVs,  and  fi 


This  chapter  briefly  reviews  the  concepts  of  Linear  Fractional  Transformations 
(LFT’s),  Linear  Matrix  Inequalities  (LMFs),  and  the  Structured  Singular  Value,  a 
matrix  function  denoted  by  It  is  intended  only  as  an  introduction  to  familiarize 
the  reader  with  concepts  used  throughout  this  thesis,  and  is  by  no  means  a  com¬ 
prehensive  review  of  the  subjects.  For  a  more  thorough  coverage,  the  reader  should 
consult  [BDG+91],  [ZPD82],  [PZPB91],  [ZDG96],  [BEFB94],  [GNLG92],  and  the 
papers  referenced  therein. 


2.1  Linear  Fractional  Transformations 


LFT’s  are  used  extensively  in  modern  control  theory  to  describe  the  relation¬ 
ship  between  components  of  a  linear  model;  most  commonly,  of  a  plant  and  its 
controller,  or  of  a  system  and  its  set  of  possible  uncertainties.  Speciflcally,  consider 
a  matrix  M  such  that 

, ,  Mil  M\2 

M  = 

M21  M22 

Now  suppose  there  is  an  appropriately  dimensioned  block  structure  A;  (such  that 
M22A;  is  square)  that  is  related  to  M  as  shown  in  Figure  2.1. 


Figure  2.1  Lower  LFT 
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The  feedback  equations  for  this  system  are  then 

ei  =  Miidi  +  Mi2U\ 

Hi  =  M2idi  +  M22U1  (2-2) 

ui  =  Aiyi 

Note  that  this  is  very  general  and  could  represent  a  state  feedback  or  output  feedback 
system,  or  the  relationship  between  a  nominal  M  and  an  uncertainty  A/,  or  any  other 
similar  relationship.  Solving  the  above  equations  for  the  transfer  function  relating 
the  input  di  to  the  output  ei  gives  the  following 


F;(M,  Ai)  =  Mu  +  Mi2 A;(/  -  M22A,)-^M21  (2.3) 

Fi(M,Ai)  is  called  the  lower  LFT on  M  by  A;  because,  as  illustrated  in  Figure  2.1, 
the  “lower”  loop  of  M  is  the  one  closed  by  A;.  In  addition,  the  system  is  considered 
well-posed  (i.e.  it  has  a  unique  solution)  only  if  the  inverse  of  (/  —  M22A1)  exists; 
otherwise,  we  have  a  singular  problem  with  either  no  solution  or  an  infinite  number 
of  solutions. 

Similarly,  Figure  2.2  illustrates  the  relationship  leading  to  the  upper  LFT  for¬ 
mula 

Fu{M,  A„)  =  M22  +  M21  A„(/  -  Mil A„)-'Mi2  (2.4) 
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The  upper  LFT  and  lower  LFT  expressions  are  each  often  referred  to  as  simply 
LFT’s,  and  the  context  or  a  diagram  indicates  the  actual  expression  to  be  used. 

Several  interpretations  of  LFT’s  will  be  used  extensively  in  this  thesis.  The 
most  common  interpretation  and  use  of  these  expressions  is  that  the  LFT’s  are  simply 
the  closed-loop  transfer  functions  from  d\  to  ti  in  the  case  of  the  lower  LFT,  and 
from  6,2  to  62  in  the  case  of  the  upper  LFT.  In  other  words, 


=F,(M,A0 


(2.5) 


=  (2.6) 

where  M  might  represent  the  controlled  plant  and  A  could  be  the  model  uncertainties 
or  the  controller. 

Another  useful  interpretation  of  an  LFT  (for  example,  Fi{M,Ai))  is  that  it 
includes  a  nominal  Mu  independent  of  A/  but  perturbed  by  A;,  while  M12,  M21,  M22 
reflect  a  priori  knowledge  of  how  the  perturbation  will  affect  the  nominal  Mu-  A 
similar  interpretation  can  be  made  with  jP„(M,  A^),  with  M22  then  considered  as  the 
nominal  mapping. 

The  LFT  expressions  can  also  be  used  to  determine  the  transfer  function  for 
the  state  space  realization  of  a  system.  A  system  given  by 


X  =  Ax  -|-  Bu 
y  =  Cx  -t-  Du 


(2.7) 


has  a  transfer  function  G(s)  =  D  +  C{sl  —  A)  ^B.  Not  surprisingly,  this  can  also 
be  expressed  as  an  LFT 


G(5)  =  F,{ 


A  B 
C  D 


(2.8) 
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where  A  =  |/  is  now  a  matrix  representing  the  frequency  structure  of  the  state 
space  realization.  In  this  case,  A  is  more  of  a  mathematical  device  than  a  “physical” 
disturbance  or  parameter  matrix. 

Such  transfer  functions  G(s)  representing  a  system  state  space  realization  will 
often  be  denoted  as 

G  = 


It  is  important  to  note  that 


is  commonly  referred  to  as  a  coefficient  matrix, 


whereas 


is  referred  to  as  a  transfer  matrix.  Unlike  the  coefficient  matrix, 


the  transfer  matrix  implies  that  an  upper  LFT  transformation  with  A  =  ^/  has  been 
performed  on  the  coefficient  matrix,  so  that  the  only  actual  input/output  channels 
remaining  are  those  corresponding  to  the  C,  D  rows.  The  A,  B  rows  now  represent 


internal  states  of  the  transfer  matrix. 


An  important  fundamental  property  of  LFT’s  is  that  linear  interconnections  of 
LFT’s  can  be  always  be  grouped  together  into  one  combined  LFT;  that  is,  all  of  the 
uncertainties  can  be  isolated  into  one  block  diagonal  structure  affecting  a  combined 
known  fixed  system.  For  example,  consider  the  interconnection  shown  in  Figure  2.3 
consisting  of  three  components,  each  of  which  is  an  LFT. 

By  collecting  all  of  the  known  systems  together  and  collecting  all  of  the  uncer¬ 
tainties  together,  the  system  can  be  redrawn  as  in  Figure  2.4  which  is  now  an  LFT 
of  a  general  known  component  perturbed  by  a  block  diagonal  uncertainty  structure. 
This  uncertainty  model  is  commonly  referred  to  as  structured  uncertainty  and  is 
due  to  the  fact  that  the  uncertainty  of  each  component  is  independent  of  the  other 
component  input/output  channels. 

There  are  many  excellent  examples  of  the  uses  of  LFT’s  in  control  or  mathe¬ 
matics  given  in  [ZDG96],[ZPD82].  One  particular  use  relevant  to  this  thesis  concerns 
expressing  a  nonlinear  polynomial  function  of  one  or  more  indeterminate  variables 
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in  an  LFT  form,  thereby  making  it  a  linear  model  of  the  function.  For  example, 
consider  the  input /output  relation 

z{0i,  O2)  =  (a  +  102  +  c0i6l)w  Xw  (2.10) 


where  a,  6,  and  c  are  constants,  and  9i  and  62  are  indeterminate.  In  order  to  express 
X  as  an  LFT  in  terms  of  and  62,  a  diagram  showing  the  input /output  relations 
with  each  9  should  first  be  drawn,  denoting  the  inputs  and  outputs  of  the  0’s  as  ?/’s 
and  u’s.  For  our  example,  this  is  shown  in  Figure  2.5. 


Figure  2.5  Block  diagram  of  X  showing  input/output  relations 


The  outputs  y  and  z  can  now  be  written  in  terms  of  the  inputs  u  and  w, 
without  referring  to  the  indeterminate  0’s.  This  essentially  pulls  out  the  0’s  into  one 
indeterminate  block  diagonal  structure.  Thus,  for  this  example: 


M 


/* 

- 

N 

yi 

0 

0 

0 

1 

Ui 

2/2 

1 

0 

0 

0 

U2 

2/3 

0 

c 

0 

b 

U3 

.2; 

0 

0 

1 

a 

w 

(2.11) 


2-6 


Therefore, 


2r  :=  Gw  =  -F«(M,  0)rti,  where  0  = 


01  0 
0  02^2 


where  /„  represents  an  n  x  n  identity  matrix. 


(2.12) 


This  technique  will  be  used  in  Chapter  IV  to  rewrite  several  polynomial  in¬ 
put/output  expressions  as  LFT’s.  These  individual  LFT’s  are  interconnected  and 
will  then  be  gathered  into  one  combined  LFT  perturbed  by  a  structured  uncertainty 
block.  As  we  will  see,  this  enables  us  to  model  many  time- varying  systems,  including 
the  aircraft  under  study. 


2.2  Linear  Matrix  Inequalities 

A  matrix  inequality  is  any  constraint  of  the  form 

A{x)  <  0 


(2.13) 


where 

•  X  —  (xi,  •  •  • ,  a:/^)  is  a  vector  of  free  design  variables,  and 

•  A{x)  is  a  symmetric  matrix. 

Note  that  “<  0”  indicates  that  A  must  be  negative  definite,  and  thus  its 
eigenvalues  are  all  negative.  Similarly,  an  expression  such  as  V  >  0  would  indicate 
that  X  is  positive  semidefinite  so  that  all  its  eigenvalues  are  positive  or  equal  to  zero. 
This  convention  will  be  used  throughout  this  thesis.  Note  also  that  the  form  of  the 
inequality  in  (2.13)  is  really  quite  general,  since  equations  can  be  manipulated  and, 
if  required,  augmented  with  dummy  variables  to  end  up  in  this  form. 

Linear  matrix  inequalities  (LMI’s)  are  a  special  class  of  matrix  inequalities 
where  A{x)  depends  affinely  on  the  design  variables  x.  That  is 

A.(a;)  :=  Aq  -f-  x^A-i  T  ■  ■  ■  F  ^NAj^f  <  0  (2.14) 
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where  Ao,  Ai, . . . ,  Ajv  are  all  given  symmetric  matrices. 

A  system  of  LMI’s  is  then  a  set  of  M  matrix  constraints  given  by 

Ai{x)  <  0,  where  ^  =  1, . . . ,  M  (2.15) 

This  system  can  be  treated  as  a  single  LMI  by  replacing  it  with 

A{x)  =  diag(Ai(a;), . . . ,  Am{x))  <  0  (2.16) 

LMI’s  are  numerically  appealing  since  finding  a  solution  x  to  satisfy  the  above  equa¬ 
tion  is  a  convex  feasibility  problem  for  which  efficient  interior-point  algorithms  are 
now  available  [NN94],[BEFB94],[GNLC92].  As  explained  in  these  references,  such 
algorithms  first  define  an  objective  function  which  is  finite  within  the  feasible  set.  In 
this  way,  the  feasibility  problem  is  transformed  into  a  convex  optimization  problem. 
The  solution  chosen,  if  any  exists,  then  corresponds  to  the  minimum  value  of  the 
objective  function,  defined  as  the  analytic  center  of  the  LMI  [BEFB94]. 

A  basic  example  of  an  LMI  is  the  Lyapunov  inequality 

A^X  -h  A:A  +  Q  <  0,  where  X  =  (2.17) 

In  such  cases,  the  variables  x  consist  of  the  free  entries  of  the  unknown  sym¬ 

metric  matrix  X. 

Many  theorems  and  expressions  using  LMI’s  are  derived  by  making  use  of  the 
Schur  complement  to  transform  equations  into  negative  definite  2x2  block  matrices 
which  are  then  LMI’s.  Specifically,  consider  the  block  matrix 


Then  T  <  0  if  and  only  if 


5  <  0 

P  -  <  0 


(2.19) 


In  this  case,  P  —  MS~^M^  is  referred  to  as  the  Schur  complement  of  S.  Note  that 
P,  S  and  M  can  themselves  be  matrices  or  matrix  expressions.  For  example,  consider 
the  Riccati  inequality 

A^X  +  +  XBB^X  +  Q  <  0,  where  X  =  X^  G  (2.20) 

Using  Schur  complements,  this  can  be  rewritten  as  the  LMI 

A^X  +  XA  +  Q  XB 

<0  (2.21) 

B^X  -I 

where  X  is  again  an  unknown  symmetric  matrix  consisting  of  the  variables  x. 

In  fact,  since  5  <  0  in  (2.19),  it  can  also  represent  a  matrix  formed  by  use  of 
the  Schur  complement.  An  important  example  of  this  is  the  Bounded  Real  Lemma 
(BRL)  for  continuous-time  systems,  which  is  used  to  turn  Hex,  constraints  into  LMI’s. 

Continuous-Time  BRL:  Consider  a  transfer  function  T(s)  =  D -\-C (si — A)~^ B . 
The  following  statements  are  then  equivalent: 

1.  A  is  stable  and  ||r(5)||oo  <  7- 

2.  There  exists  a  solution  X  >  0  to  the  following  inequalities 
a(D)  =  ||T)||oo  <  7 

A^X  +  XA  +  +  j{XB  +  -  D^D)-\B^X  +  ^D^C)  <  0 

(2.22) 
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3.  There  exists  a  solution  X  >  0  to  the  LMI 

A^X  +  XA  XB 

B^X  -7/  <  0.  (2.23) 

C  D  -7/ 

One  can  clearly  appreciate  the  power  of  LMI’s  here,  noting  that  the  two  conditions 
of  item  2,  including  the  rather  complex  Algebraic  Riccati  Inequality  (ARI),  can  be 
expressed,  after  some  manipulations,  as  the  LMI  of  item  3;  also,  the  LMI  version 
can  be  solved  efficiently  since  it  is  a  convex  programming  problem. 

This  BRL  forms  the  basis  of  the  proofs  for  LMI-based  Hqo  control,  and  thus 
for  the  theorems  and  results  to  be  presented  in  Chapter  III. 

In  general,  there  are  three  types  of  LMI  problems: 

1.  determining  a  feasible  solution  x  (or  X)  for  some  LMI  constraints  given  by 
L{x)  <  0. 

2.  minimization  of  a  convex  objective  given  some  LMI  constraints;  that  is, 

minimize  c^a;  subject  to  L{x)  <  0.  (2.24) 


3.  solving  the  generalized  eigenvalue  minimization  problem  (GEVP)  given  some 
LMI  constraints;  that  is 


A{x)  <  \B{x] 

minimize  A  subject  to  <  B{x)  >  0 

C{x)  <  0. 


(2.25) 


The  first  two  types  of  problems  are  convex,  while  the  GEVP  problem  is  quasi-convex. 
The  feasibility  and  GEVP  problems  will  be  referred  to  in  Chapter  III  and  used  to 
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obtain  the  gain-scheduled  controller.  All  of  these  three  types  of  problems  can  now 
be  addressed  using  the  new  LMI  Control  Toolbox  in  MATLAB  [GNLC92],[MAT]. 

As  a  final  example  of  the  power  of  LMFs,  recall  that  in  the  standard  Riccati- 
based  [DGKF89]  method  of  computing  controllers,  the  Algebraic  Riccati  Equa¬ 
tion  (ARE) 

A^X  +  XA  +  X(j-^BiB[  -  B2B^)X  +  CfCi  =  0  (2.26) 

is  used  to  find  a  stabilizing  solution  X^o,  which  is  then  used  to  determine  the  “cen¬ 
tral”  Hoo  controller  (note  that  D  =  0  is  assumed  here  for  simplicity).  All  suitable 
controllers  can  then  be  parametrized  via  LFT’s  built  around  the  central  controller 
using  a  free  parameter  Q  [DGKF89].  A  problem  with  this  approach  is  that  the  plant 
must  be  regular,  and  the  resulting  controllers  are  often  of  high  order,  especially 
when  using  Q-parametrization,  since  there  is  no  clear  connection  between  Q  and  the 
controller  properties. 

In  the  LMI-based  approach,  the  ARE’s  are  replaced  by  ARFs,  and  the  solution 
set  of  all  these  inequalities  is  used  to  parametrize  the  Hoo  controllers.  As  a  result, 
the  plant  need  not  be  regular,  and  the  resulting  controller  orders  and  properties  can 
be  shown  to  depend  on  two  matrix  parameters  R  and  S  [GA94].  In  addition,  the 
LMI  constraints  can  be  altered  or  other  LMI  constraints  can  be  added  to  further 
define  or  modify  the  problem.  Since  a  set  of  LMI’s  remains  an  LMI,  this  is  an 
elegant  way  to  ensure  the  problem  remains  convex  and  optimizable.  This  approach 
is  taken  in  Chapter  III  to  formulate  the  problem  in  order  to  be  able  to  solve  for  the 
gain-scheduled  controller. 

2.3  The  Structured  Singular  Value 

Unstructured  uncertainty  is  generally  modelled  as  some  norm-bounded  A  € 
Hoo  perturbing  a  nominal  plant  or  system,  as  shown  in  Figure  2.6. 
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Figure  2.6  Modelled  uncertainty  A  in  a  system 


The  problem,  then,  becomes  one  of  either  finding  the  largest  A  the  system  can 
handle  or  of  finding  the  best  controller  K  to  maximize  the  allowable  A  and,  in  so 
doing,  maximize  the  stability  margins.  This  is  the  robust  stability  problem,  and  is 
addressed  by  using  the  Small  Gain  Theorem  to  derive  tests  for  robust  stability. 

However,  this  does  not  directly  address  robust  performance,  in  which  we  would 
like  to  maintain  both  robust  stability  and  nominal  performance  throughout  the  en¬ 
velope  modelled  by  the  uncertainty  A.  In  other  words,  we  want  the  plant  not  only 
to  be  stable  to  some  uncertainties,  but  also  to  perform  to  a  desired  level  even  when 
perturbed.  There  are  several  possible  ways  to  do  this,  and  these  all  involve  maximiz¬ 
ing  both  a  test  for  robust  stability  and  another  test  for  performance,  and  performing 
some  tradeoffs  for  the  optimal  mix.  Combining  both  subproblems  as  one  test  im¬ 
plies  a  structure;  for  strictly  Hoo  problems,  this  has  led  to  the  development  of  the 
Structured  Singular  Value,  called  SSV  or  p,  in  order  to  optimize  this  combined  test. 

Looking  again  at  Figure  2.6,  we  see  that  the  transfer  function  from  6,2  to  62  is 
given  by 

=  F„(M,Ai)  (2.27) 

where  M  is  as  defined  in  (2.1).  This,  of  course,  assumes  that  the  LFT  is  well-posed 
and  that  Ai  has  a  block  structure  compatible  with  Mu.  In  this  case,  M22  is  the 
nominal  map  and  the  rest  of  M  reflects  how  the  norm-bounded  perturbation  Ai 
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affects  M22-  Now  consider  a  structure 


A  = 


Ai  0 
0  A2 


(2.28) 


As  mentioned  in  the  last  section,  such  an  uncertainty  block  is  considered  a  structured 
uncertainty  since  it  results  from  combining  all  the  uncertainties  at  different  points  of 
a  system  into  one  block.  In  this  case,  since  we  are  interested  in  robust  performance, 
A2  could  represent  a  fictitious  uncertainty  associated  with  the  performance  channel; 
that  is,  from  62  to  ^2.  The  resulting  interconnection  is  shown  in  Figure  2.7.  In 
fact,  if  we  have  multiple  uncertainties  in  our  system  (or  even  multiple  performance 
objectives),  our  uncertainty  block  structure  could  be  composed  of  several  more  A’s 
linked  to  the  appropriate  input-output  channels. 


In  general,  there  are  two  types  of  uncertainty  blocks:  repeated  scalar  blocks, 
and  full  blocks.  For  S  repeated  scalar  blocks  6i,  each  of  dimension  rs  X  rs,  and  F 
full  blocks  Aj,  we  can  define  a  general  structured  uncertainty  A  G  as 


(2.29) 


2-13 


For  dimensional  consistency,  if  M  G  then  ri  +  Ylf=i  '^j  —  The  norm- 

bounded  subset  of  A  can  be  defined  as 


=  {A  G  A  :  a(A)  <  1} 


(2.30) 


where  cr(A)  is  defined  as  the  maximum  singular  value  of  A. 

The  structured  singular  value  /x  is  a  function  created  to  get  a  measure  of  the 
effect  of  the  smallest  perturbation  A  G  A  for  which  a  system  becomes  singular  and 
therefore  unstable.  It  is  defined  as 

min  {d(A)  :  A  G  A,  det(7  -  MA)  =  0} 

unless  I  —  MA  is  always  nonsingular,  in  which  case  /ia  :=  0. 

For  example,  in  our  earlier  problem  of  (2.27)-(2.28)  where  A  is  a  structured 
uncertainty  comprised  of  an  actual  uncertainty  and  a  fictitious  performance  “uncer¬ 
tainty”,  we  can  observe  that  minimizing  fi  will  maximize  the  size  of  the  allowable 
A,  which  is  exactly  what  we  would  like  to  achieve!  In  fact,  this  leads  us  to  the  main 
theorem  for  the  use  of  fi  in  linear  system  robustness  analysis; 

Main  Loop  Theorem:  The  following  equations  are  equivalent: 


)«A  W  <  1  1 


MAi(Mii)  <  1,  and 
AiG^a,  Ma.(WAi))<1 


(2.32) 


Proof.  See  [ZDG96]. 

The  first  equation  on  the  right-hand  side  is  p  of  Mu  with  respect  to  the  uncer¬ 
tainty  Ai.  This  is  basically  equivalent  to  the  Small  Gain  Theorem,  and  guarantees 
well-posedness  and  robust  stability  if  less  than  1.  The  second  equation  on  the  right- 
hand  side  is  fi  of  F)(M,  Ai)  with  respect  to  the  fictitious  uncertainty  A2.  This  is 
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basically  a  measure  of  the  inverse  of  the  performance  of  the  perturbed  system.  The 
point  is  that  if  /ua  is  made  less  than  1,  we  can  then  guarantee  robust  stability  and 
some  level  of  performance;  in  other  words,  robust  performance. 

Expression  (2.31)  for  is  very  hard  to  evaluate  in  practice;  however,  for  the 
extreme  case  where  A  G  then 

fx^{M)  =  a{M)  (2.33) 

This,  then,  is  an  overbound  to  fj,  for  less  extreme  cases  where  A  is  a  mixture  of  full 
blocks  and  repeated  scalar  blocks;  i.e., 

Ma(M)  <  a{M)  (2.34) 

In  order  to  close  the  gap  in  the  overbound,  positive  definite  similarity  transformations 
on  M  can  be  found  that  affect  the  value  of  a  but  do  not  affect  the  value  of  /u.  This 
is  given  by  the  following  theorem 

Theorem  2.3.1:  For  all  Z)  G  D  and  A  G  A,  where  DA  =  AD  and 

D  =  {diag[Di, . . . ,  Ds,  dilm^ ,  • . . ,  dplmp]  :  Di  G  =  D*  >  0,  dj  G  dj  >  0} 

(2.35) 

then 

IIa{M)  =  ha{DMD-^)  (2.36) 


Proof.  See  [ZDG96]. 

It  is  important  to  note  that  the  transformations  or  “£)-scales”  are,  like  the 
corresponding  uncertainties,  of  two  types.  For  repeated  scalar  uncertainties  the  D- 
scales  are  full  blocks;  whereas  for  full  block  uncertainties  the  D-sc&les  are  scalars. 
Thus,  if  we  would  want  to  normalize  the  scaling  for  an  uncertainty  A^  by  the  scaling 
for  an  uncertainty  A;,,  it  becomes  necessary  to  treat  the  uncertainty  Aj  as  a  full 
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block  in  order  for  its  corresponding  Z)-scales  to  be  scalar.  This  is  done  in  Chapter 
III  in  the  D-K-D  scheme  presented  there,  even  though  in  that  case  the  uncertain¬ 
ties  were  known  to  be  repeated  scalar  blocks.  This  unfortunately  introduces  some 
conservativeness,  since  we  are  then  allowing  for  uncertainties  that  are  known  not  to 
exist. 

Returning  to  our  discussion,  combining  the  results  of  Theorem  2.3.1  and  (2.34), 
it  follows  that 

/<a(V)  < =  d'zD  HDMD-')  (2.37) 

At  this  point,  it  should  be  pointed  out  that  n  is  not  necessarily  constant.  In 
the  previous  development,  we  assumed  that  both  M  and  A  were  constant,  and  for 
such  cases,  would  be  constant;  however,  as  is  often  the  case,  M  or  A  (or  both)  is 
usually  a  system  matrix  varying  with  frequency.  For  such  cases,  n  will  be  a  frequency 
dependent  function.  The  previous  development  remains  valid,  except  that  in  order 
to  determine  the  maximum  value  for  jx,  we  must  then  evaluate  it  over  all  frequencies 
for  its  supremum.  To  handle  these  more  general  cases,  (2.37)  must  then  be  modified 
slightly  as  follows: 

a{DMD-^)  (2.38) 

Since  we  don’t  really  know  how  to  calculate  /x,  the  upper  bound  p,  is  what  is 
used  to  approximate  it.  Determining  the  D-scales  is  a  convex  programming  problem, 
so  these  can  be  solved  for  efficiently  [PD93],[BDG'*'91].  The  value  of  fx  can  then 
be  approximated  over  all  frequencies  to  determine  its  frequency  response  and  its 
maximum.  There  is  also  a  lower  bound,  but  since  it  is  nearly  as  difficult  to  evaluate 
as  IX,  it  is  not  commonly  used.  In  fact,  for  certain  uncertainty  block  structures, 
p  has  been  shown  to  be  equal  to  fi  [PD93].  Most  importantly,  if  the  uncertainty 
is  composed  of  1  to  3  full  blocks  (and  no  repeated  scalars)  then  the  upper  bound 
is  equal  to  fx.  It  has  also  been  shown  to  remain  relatively  tight  for  4  or  more  full 
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blocks.  Thus,  using  a  program  to  solve  for  the  D-scales  such  as  the  /x-Toolbox  in 
MATLAB  [BDG‘''91],[MAT],  we  can  approximate  both  the  frequency  response  and 
the  maximum  value  of  fi  for  a  system.  This  is  called  /^-analysis,  and  it  can  be  used 
to  guarantee  that  a  general  system  exhibits  robust  performance  by  demonstrating 
that  the  peak  value  of  fi  <  1. 

However,  a  more  common  problem  is  that  of  finding  a  stabilizing  controller 
that  would  maximize  robust  performance  and  thus  minimize  pL.  This  is  called  //- 
synthesis.  For  a  general  system  like  the  one  shown  in  Figure  2.8,  we  would  therefore 
like  to  find  a  stabilizing  K  such  that  the  peak  value  of  jw  <  1.  Noting  that  the 
transfer  function  from  the  inputs  d  to  the  outputs  e  is  Ted  =  we  can  then 

express  the  //-synthesis  problem  as 

=  “Sf  S(DT^D-^) 

(2.39) 

Therefore,  a  tight  approximation  to  our  desired  //-synthesis  robust  performance  prob¬ 
lem  is  given  by 

ofo  \\DT,.D-'\U  (2.40) 

Given  rational  functions  for  the  D-scales,  this  expression  is  now  a  standard,  albeit 
scaled,  H^o  problem  which  can  be  solved  using  Riccati-based  or  LMI-based  methods 
for  the  optimal  controller.  This  scaled  problem  is  illustrated  in  Figure  2.9. 

However,  solving  for  both  the  D-scsles  and  K  is  not  a  jointly  convex  prob¬ 
lem.  The  most  popular  //-synthesis  approach  is  the  so-called  D-K  iteration  method, 
summarized  in  the  following  steps; 

1.  Set  D  =  I  initially. 
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2.  Solve  the  Hoo  optimization  problem  \\DTedD  ^||oo  for  the  minimizing,  stabiliz¬ 
ing  controller  K.  Note  that  Ted  =  Fi{PiK). 

3.  Use  this  K  and  determine  the  D-scales  to  minimize  Ji  given  by 

n'to  \\DT^D-^\U.  (2.41) 

4.  Fit  the  D-scales  to  rational  transfer  functions  to  closely  approximate  ji. 

5.  Repeat  steps  2  to  4  until  ji  no  longer  decreases. 

Note  that  this  procedure  is  not  guaranteed  to  converge,  but  since  it  works  well 
in  practice  it  is  implemented  in  the  /i-Toolbox  in  MATLAB.  Also,  the  discussion 
of  dealt  specifically  with  complex  uncertainties.  For  real  or  mixed  uncertainties, 
there  are  similar  theoretical  results  [BDG+Ql];  unfortunately,  while  )W-analysis  has 
been  implemented  for  both  complex  and  real  uncertainties,  /x-synthesis  has  not.  This 
is  unfortunate,  since  it  forces  practitioners  to  assume  complex  uncertainties  for  fx- 
synthesis  even  when  some  or  all  of  these  are  known  to  be  real,  resulting  in  a  much 
more  conservative  design  than  would  otherwise  be  necessary.  This,  for  example,  is 
the  case  for  the  aircraft  controller  designed  in  Chapter  IV. 
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III.  Gain  Scheduling  Theory 

This  chapter  presents  the  background  theory  of  gain  scheduling  for  linear  pa¬ 
rameter  varying  systems  as  developed  principally  by  Packard  [Pac94],  and  Apkarian 
and  Cabinet  [AG95].  Most  of  the  following  development  is  taken  with  only  nota- 
tional  and  other  minor  modifications  from  [AG95]  and  [Pac94].  While  the  following 
development  is  given  only  for  continuous-time,  the  cited  references  also  detail  the 
discrete-time  approach. 

3.1  Linear  Parameter- Varying  Systems 

Much  of  Hoo  theory  to  this  point  has  involved  the  synthesis  of  controllers  for 
linear  time-invariant  (LTI)  plants.  This  applies  well  for  a  host  of  problems  where 
the  plants  remain  relatively  constant  over  time.  However,  in  many  applications,  the 
plant  varies  significantly  over  time.  Such  linear  time-varying  (LTV)  plants  can  be 
expressed  with  state-space  equations  of  the  form 

x{t)  =  A{t)x{t) B{t)u{t) 
y{t)  =  C{t)x{t)  -\-  D{t)u{t) 

where  not  only  the  states  x(t)  and  the  inputs  u{t)  vary  with  time,  but  also  the  state- 
space  matrices  A(f),  B{t),  C'(f),  and  D[t).  Further,  a  large  class  of  these  LTV  plants 
can  be  expressed  as  systems  of  the  form 

x(t)  =  A{6{t))x(t)  -\-  B{0{t))u{t) 
y{t)  =  C{6{t))x{t)  -I-  D{B{t))u{t) 

where  9{t)  is  a  vector  of  time- varying  parameters  and  A(0),  B{d).,  C{9)^  D{9)  are 
now  known  fixed  functions  of  0{t).  The  difference  lies  in  the  fact  that,  while  the 
state-space  matrices  of  (3.1)  are  known  functions  of  time,  the  state-space  matrices  of 
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(3.2)  are  known  functions  of  0(t).  0[t)  must  be  measured  real-time  and  is  not  known 
a  priori,  although  its  range  of  parameter  values  is  often  known. 

Such  plants  are  now  commonly  referred  to  as  linear  parameter-varying  (LPV) 
systems,  and  can  describe  many  types  of  systems  such  as  aircraft,  robotics,  missiles, 
etc.  In  the  case  of  an  aircraft,  its  stability  derivatives  and  thus  its  state-space  model 
depends  on  parameters  such  as  Mach  number,  altitude,  angle  of  attack,  or  dynamic 
pressure;  these  parameters  in  turn  change  over  time. 

In  many  cases,  A{6),  B{6),  C{6),  D{9)  can  be  expressed  as  a  linear  fractional 
function  of  9{t).  In  other  words,  we  can  express  A{6),  B(6),  C{6),  and  D{6)  as  an 
LFT  relating  6{t)  to  an  LTI  plant  P  consisting  of  “nominal”  values  A,  B,  C,  and  D. 
Thus,  using  standard  state-space  notation, 


P(s) 


A 

B 

C 

D 

(3.3) 


The  time-varying  parameter  vector  can  be  represented  as 

e  =  {9r,---,9K)^^^  (3.4) 

Referring  to  Figure  3.1,  0  is  a  block  diagonal  time- varying  operator  specifying  how 
9  enters  the  plant  dynamics.  Specifically, 


0  =  blockdiag(^i/ri ,  •  •  • ,  9kItk)  (3-5) 

where  rj  >  1  whenever  the  parameter  9i  is  repeated  [Doy85].  The  set  of  all  0  is 
given  by 

A  :=  {0  :  9i{t)  G  (3.6) 

Note  that  A  is  traditionally  referred  to  as  an  uncertainty  structure.  In  fact,  if  we 
also  wanted  to  model  some  uncertainty  in  the  plant,  the  structure  of  A  could  be 
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modified  to  account  for  this.  This  will  be  expanded  upon  later.  Note  also  that 
0  e  3?  since  the  actual  plant  variations  are  real,  not  complex;  however,  since  we 
will  eventually  be  applying  the  Small  Gain  Theorem  to  solve  for  our  controller,  the 
solution  will,  in  fact,  allow  for  complex  variations  in  the  plant.  This  results  in  a 
much  more  conservative  solution  and  less  performance  than  could  likely  be  achieved. 


Figure  3.1  LPV  Plant 


The  feedback  equations  for  the  plant  LFT  can  now  be  written  as 


P(s) 


Thus,  for  a  given  time  t,  the  LPV  plant  defines  an  LTI  plant  whose  transfer  function 
is  given  by  the  upper  LFT 


=  t;(p,0) 

.  y . 


(3.9) 
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or,  explicitly, 


e  ^  Ped  Pen  PeS  ©[/ _  (3.10) 

y  Pyd  Pyu  Py6  .  ^ 

Using  the  LPV  model  of  the  plant,  traditional  techniques  based  on  the  Small 
Gain  Theorem  could  now  be  used  to  find  a  single  robust  LTI  controller  for  the  vary¬ 
ing  plant  [Doy85].  Unfortunately,  if  the  plant  undergoes  large  deviations,  this  will 
generally  be  very  conservative  and  exhibit  poor  performance.  In  fact,  a  stabilizing 
controller  may  not  even  be  feasible.  A  better  approach  is  to  design  a  parameter- 
dependent  controller,  which  will  then  allow  it  to  vary  as  the  plant  varies.  This  is 
illustrated  in  Figure  3.2. 


Figure  3.2  LPV  control  structure 


The  feedback  equation  for  the  gain-scheduled  LPV  controller  is  then 

u  =  Fi{K,Q)y  (3.11) 

where  Fi{K,Q)  is  the  controller  at  a  given  time  t,  and  K  specifies  the  relationship 
between  the  changing  parameters  and  the  gain- scheduled  controller.  Since  0  is  based 
on  the  known  parameter  variations  of  the  plant,  the  only  unknown  is  the  LTI  K  given 
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by 


K{s)  = 


Kuy 

Key  Ko0 


(3.12) 


de  =  Qee 


(3.13) 


Thus,  once  K  has  been  found,  the  LPV  gain-scheduled  controller  could  then  be 
physically  implenaented.  The  controller  would  then  be  continually  updated  with  in¬ 
flight  measurements  of  6,  and  (3.11)  gives  the  rule  to  produce  the  updated  control 
input  u. 

Finally,  note  that  the  closed-loop  transfer  function  from  the  disturbance  d  to 
the  controlled  output  e  is  given  by 


^  =  nd{P,K,Q)  =  Fi{F4P,e),Fi{K,Q)).  (3.14) 

a 

Note  also  that,  without  any  0,  we  recover  the  transfer  function  associated  with  LTI 
systems;  namely  Ted  =  Fi{Pyd(). 


3.2  Hoo  Control  of  LPV  Systems 

Given  an  LTI  plant  P{s),  the  usual  Hqo  problem  requires  finding  an  internally 
stabilizing  LTI  controller  K(s)  such  that 

||f’,(-P.^)IU<7  (3.15) 

where  7  >  0  is  some  desired  performance  level.  However,  in  the  gain-scheduled 
version  of  this  problem,  both  the  plant  and  the  controller  are  LPV  vice  LTI.  The 
objective  becomes  one  of  robust  performance;  in  other  words,  to  guarantee  some 
closed- loop  performance  7  for  all  admissible  parameters  6.  The  Hoo  control  problem 
is  now  required  to  find  a  controller  K{s)  such  that  the  LPV  controller  Fi{K,  0)  meets 
the  following  conditions 


3-5 


1.  the  closed- loop  system  must  be  internally  stable  for  all  admissible  parameters 

6  such  that 

<  1  (3.16) 

2.  the  closed-loop  system  satisfies 

\\Ted{P,K,e)\\oc  <lf  (3.17) 

l|0||oo  <  1/7 

Equation  (3.16)  comes  from  the  Small  Gain  Theorem  and  restricts  the  parameter 
range  of  0  to  a  ball  of  radius  I/7.  This  implies  no  loss  of  generality  since  the  param¬ 
eters  9  can  simply  be  scaled  to  comply  with  this  requirement.  For  example,  in  the 
implementation  described  in  Chapter  IV,  the  parameter  variations  have  been  nor¬ 
malized  such  that  7  <  1  would  imply  that  the  above  two  conditions  have  been  met. 
Following  the  convention  of  [AG95],  solutions  to  the  aforementioned  ifoo  problem 
for  LPV  systems  will  be  called  7-suboptimal  gain-scheduled  controllers. 

As  shown  in  Figure  3.2,  the  gain-scheduling  problem  as  stated  has  a  time- 
varying  parameter  block  0  entering  both  the  plant  and  the  controller.  In  order  to 
make  use  of  the  Small  Gain  Theorem,  the  two  parameter  blocks  can  be  combined 
into  one  single  uncertainty  block,  as  illustrated  in  Figure  3.3,  resulting  in  a  new  plant 
Pa. 

The  feedback  equations  for  Pa  can  be  written  as 


PM 

de 

ee 

0 

0 

Ir 

de 

e 

— 

0 

P{s) 

0 

d 

y 

Ir 

0 

0 

u 

y 

u 
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Figure  3.3  Modified  LPV  control  structure 
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As  a  result,  the  closed-loop  transfer  function  between  the  exogenous  input  d  and  the 
controlled  output  e  is  now 


Ted{P,  K,  0)  =  FuiFiiPaJ<),  A  ©  A),  where  A  0  A  = 


0  0 
0  0 


(3.19) 


In  state-space,  the  LTI  plant  P  can  be  written  as 


P(s)  = 


A 

Be 

Bd 

Bu 

Ce 

Dee 

Ded 

D$ii 

C-e 

Dee 

Ded 

Deu 

Cy 

Dye 

Dyd 

D  yu 

Thus,  the  augmented  LTI  plant  Pa  is 


P.M 


’  A 

0 

Be 

Bd 

Bu 

0  ’ 

0 

0 

0 

0 

0 

Ir 

Ce 

0 

Dee 

Ded 

Deu 

0 

Ce 

0 

Dee 

Ded 

Den 

0 

Cy 

0 

Dye 

Dyd 

Dyu 

0 

0 

Ir 

0 

0 

0 

0 

Similarly, 


Ak 

BKy 

Bxe 

K{3)  = 

Cku 

Dxuy 

Dxue 

Ckb 

DKey 

Dxee 

(3.20) 


(3.21) 


(3.22) 


For  a  problem  with  n  states,  m2  control  inputs  u,  p2  measured  outputs  y,  pi  exoge¬ 
nous  inputs  d  and  controlled  outputs  e,  time- varying  operator  0  of  dimension  r,  and 
a  controller  K  of  order  k,  the  problem  dimensions  are 


(3.23) 
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Note  that  0  and  the  transfer  function  from  the  disturbances  d  to  the  controlled 
outputs  e  have  been  considered  square.  This  simplifies  the  notation  and  can  be 
easily  met  by  augmenting  the  problem  with  rows  and/or  columns  of  zeros. 

There  are  also  three  assumptions  required  to  solve  for  the  gain-scheduled  con¬ 
troller: 

1.  {A,Bu,Cy)  must  be  stabilizable  and  detectable, 

2.  Dyu  =  0,  and 

3.  DyB  =  0  or  Dbu  =  0. 

The  first  assumption  is  standard,  and  is  necessary  and  sufficient  to  guarantee  the 
existence  of  a  stabilizing  controller.  The  second  assumption  is  not  required,  but 
greatly  simplifies  the  calculations  in  the  rest  of  this  chapter.  Appendix  A  of  [A1195] 
outlines  the  shifting  technique  that  can  be  used  to  lift  this  restriction.  The  third 
assumption  is  sufficient  but  not  necessary  to  guarantee  that  the  LPV  controller 
is  causal  and  well-posed  [AG95].  Many  problems  meet  this  condition  since  the  y 
measurement  equation  is  often  independent  of  the  parameters  0,  or  9  is  independent 
of  the  control  input  u.  This  assumption  and  the  previous  ones,  for  example,  all  hold 
for  the  problem  described  in  Chapter  IV.  Note  that  the  last  assumption  can  also  be 
removed  by  imposing  a  well-posedness  constraint  on  the  problem  [AG95]. 

We  now  have  a  more  classical-looking  robust  performance  problem,  with  a 
norm-bounded  repeated  uncertainty  block  A  0  A,  an  LTI  controller  K,  and  a  new 
LTI  plant  Pa  consisting  of  the  original  plant  P  augmented  with  the  interconnections 
u  and  y  between  K  and  0.  It  is  precisely  because  of  these  extra  interconnections 
that  this  problem  differs  from  the  more  classical  Hoo  control  problem,  since  they 
ensure  that  the  controller  changes  with  exactly  the  same  parameter  variations  as  the 
plant.  The  only  difference  lies  in  how  much  the  controller  changes  in  response  to 
the  variations;  this,  as  we  will  see  shortly,  will  be  optimized  using  LMI’s  to  find  the 
optimal  scalings  or  weights  which  will  dictate  the  magnitude  of  the  response  to  the 
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parameter  variations.  Those  optimized  scales  in  turn  will  then  be  used  to  solve  for 
the  7-suboptimal  controller  K. 

Now  that  the  problem  has  been  reformulated  to  have  a  single  uncertainty  block, 
the  Small  Gain  Theorem  can  be  used  to  determine  a  sufficient  condition  for  robust 
performance  or,  equivalently,  for  the  existence  of  a  gain-scheduled  controller  [Doy85]. 
Specifically,  this  leads  to  the  following  theorem  from  [AG95]. 

Theorem  3.2.1:  Consider  the  set  of  all  positive  definite  similarity  scalings  com¬ 
muting  with  A 


K 


La  =  {L  >  0  :  LQ  =  0L,V0  G  A}  C  9?’'^’’,  where  r  =  (3.24) 

i=l 


Given  La  ,  the  set  of  scalings  commuting  with  A  0  A  is  then 


Ta®a  =  s 


L  = 


L\  L2 
Ll  Ls 


>  0  :  Li^L^  G  LajL2Q  —  0Xr25^€)  G  A 


(3.25) 


If  there  exists  a  scaling  matrix  L  G  La®a  and  an  LTI  control  structure  K  such  that 
the  nominal  closed-loop  system  Fi{K,  0)  is  internally  stable  and  satisfies 


0 

0  I 


Fl{Pa,K) 


0 

0  I 


<  7 


00 


(3.26) 


then  Fi{K,Q)  is  a  7-suboptimal  gain-scheduled  Hoo  controller. 

Proof:  See  [AG95]. 

Note  that,  in  (3.26),  L  G  since  it  is  the  scaling  matrix  for  A  ©  A;  and 

I  G  9?*^*^^^  because  the  controlled  outputs  e  and  the  exogenous  inputs  d  are  not 
scaled.  Their  weights,  if  any,  must  be  incorporated  into  Pa.  The  problem  is  now  one 
of  calculating  L  and  K  such  that  (3.26)  holds. 
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3.3  Solving  the  General  Scaled  Hoo  Problem 

The  Bounded  Real  Lemma  (BRL)  can  be  used  to  transform  scaled  Hoo  op¬ 
timization  problems  into  a  set  of  LMI  constraints.  For  convenience,  it  is  stated 
below. 

Continuous-Time  Scaled  BRL:  Consider  an  arbitrary  uncertainty  structure  A, 
and  the  associated  scaling  set  La  as  defined  in  (3.24).  Also  consider  a  square  transfer 
function  T{s)  =  C{sl  —  A)~^B  +  D.  The  following  statements  are  then  equivalent; 


1.  A  is  stable  and  there  exists  L  G  La  such  that  \\L^^'^T{s)L  ^/^||oo  <  7- 

2.  There  exists  solutions  X  >  0  and  L  G  La  to  the  matrix  inequality 


A^X  +  XA 

XB 

CT 

B^X 

-jL 

C 

D 

-7I- 

<  0.  (3.27) 


Proof:  See  [GA94]. 

Now  consider  a  proper  LTI  plant  G{s)  where 

d 
u 


e 

'  Ged 

Geu 

y 

Gyd 

1 

0 

(3.28) 


Given  a  desired  objective  7  >  0,  an  arbitrary  uncertainty  structure  A,  and  the 
associated  scaling  set  La  as  defined  in  (3.24),  the  general  scaled  Hoo  problem  is 
one  of  finding  L  G  La  and  an  LTI  controller  K  such  that  the  closed-loop  system  is 
internally  stable  and 

\\L^I'^Fi{G,K)L-^l‘^\\oo<i.  (3.29) 


The  Scaled  BRL  can  then  be  used  to  derive  the  following  theorem  character¬ 
izing  the  solution  of  the  general  scaled  Hoo  problem. 
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Theorem  3.3.1:  Given 


A 

Bd 

Bu 

G(s)  = 

a 

Bed 

Beu 

Cy 

Byd 

0 

(3.30) 


satisfying  assumptions  1  and  2,  the  general  scaled  Hoo  problem  is  solvable  if  and 
only  if  there  exist  pairs  of  symmetric  matrices  {R,S)  G  and  {L,J)  G 

such  that 

[  AR  +  RA^  RCj  Bd 


CeR 

BJ 


-iJ  Dei 

Dli  -1L 


Nr  <  0 


(3.31) 


A^S  +  SA  SBd  Cj 
BJS  -jL  Dl 

Ce  Bed  -iJ 


Ns<Q 


R  I 
I  S 


>  0 


La,  J  ^  La,LJ  =  I 


(3.32) 


(3.33) 

(3.34) 


where  Nr  and  Ns  denote  bases  of  the  null  spaces  of  {B^,Dj^,0)  and  {Cy,Dyd,0), 
respectively. 

Further,  there  exist  suboptimal  controllers  of  order  k  if  and  only  if  (3.31)-(3.34) 
hold  for  some  [R,  S,  L,  J)  where  R,  S  satisfy  the  rank  constraint 


rank(/  —  RS)  <  k. 


(3.35) 


Proof:  See  [AG95]  and  [GA94]. 

Matrix  inequalities  (3.31)-(3.33)  are  LMI’s  in  R,  S,  L,  J  and  constraints  L  ^  La 
and  J  G  La  are  convex.  As  well,  in  the  full  order  case,  k  >  n  so  that  the  rank 
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constraint  (3.35)  is  trivially  satisfied.  However,  the  constraint  LJ  =  /  is  highly  non- 
convex,  so  that  solving  the  general  scaled  Hoo  problem  remains  a  difficult  problem. 

It  should  also  be  pointed  out  that  for  the  classical  Hoo  control  problem,  L  and 
J  can  be  set  to  I  without  loss  of  generality.  Equations  (3.31)-(3.34)  then  reduce 
to  a  set  of  three  convex  LMI’s  which  can  be  solved  numerically.  Compared  to  the 
standard  Riccati-based  method,  this  has  the  advantage  of  allowing  imaginary  axis 
invariant  zeros  and  singular  problems  arising  due  to  rank  deficiencies  in  Deu  and 
Dyd-  See  [GA94]  for  more  on  this  subject. 


3.4  Solving  the  Gain-Scheduled  Hoo  Problem 


The  gain- scheduled  Hoo  problem  can  be  viewed  as  a  specific  case  of  the  general 
scaled  Hoo  problem.  This  is  evident,  for  example,  from  comparing  (3.26)  and  (3.29). 
The  main  problem  in  solving  the  general  scaled  problem  was  the  nonconvexity  of  the 
LJ  =  I  constraint.  It  turns  out,  however,  that  due  to  the  particular  structure  of  the 
LPV  problem,  this  difficulty  can  be  completely  eliminated. 


Theorem  8.4. 1:  Given  the  LPV  system  defined  in  Sections  3.1  and  3.2  and 
satisfying  assumptions  1  and  2,  the  gain-scheduled  Hoo  problem  is  solvable  if  and 
only  if  there  exist  pairs  of  symmetric  matrices  {R,S)  G  and  (LsjJs)  G 

such  that 


Zr^N^ 


Zs  =  N^s 


AR  +  RA^ 

RCf 

BdJs 

CoR 

-ih 

DedJz 

A 

0 

(3.36) 

UL 

- 1 

1 

A^S  -h  SA 

SBd 

CJL ' 

BJS 

1 

CO 

i>Lt, 

0 

V 

(3.37) 

L3C0 

LsDed 

CO 

1 

R  I 
I  S 


>  0 


(3.38) 
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i'3  £  ^ 


L3  I 


>  0 


(3.39) 


where 

Cg 

Dgg  Dgd 

0 

CO 

Js  0 

Bd  =  [BBBd],C,= 

.  . 

?  Ded  — 

Deff  Dgd 

7^3  = 

0  / 

5  JS  — 

0  I 

(3.40) 

and  where  Nr  and  Ns  denote  bases  of  the  null  spaces  of  and 

{Cy,  Dye,  Dyd,  O),  lespectlvely. 

Further,  there  exist  suboptimal  controllers  of  order  k  if  and  only  if  (3.36)-(3.39) 
hold  for  some  (R,  S,  L3,  J3)  where  R,  S  satisfy  the  rank  constraint 


rank(7  —  RS)  <  k. 


(3.41) 


Proof:  See  [AG95]. 

Equations  (3.36)-(3.38)  are  LMI’s  in  R,S,Lz,J3  and  are  now  all  convex,  such 
that  testing  these  conditions  for  a  feasible  solution  is  a  convex  feasibility  problem. 
As  explained  in  Chapter  II,  software  using  interior-point  LMI  solvers  [NN94]  such 
as  the  LMI  Control  Toolbox  in  MATLAB  [MAT],[GNLC92]  can  then  be  used  to 
test  the  feasibility  of  the  LMI  conditions  for  some  arbitrary  7  values.  Thus,  we  can 
iteratively  solve  for  the  minimum  7  corresponding  to  the  sub-optimal  gain-scheduled 
controller. 

Alternatively,  7  can  be  solved  for  directly  since  minimizing  7  subject  to  the 
feasibility  of  these  LMI  constraints  is  a  convex  optimization  problem  [GA94].  Nr 
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and  Ns  as  defined  in  Theorem  3.4.1  can  be  written  as 


0 

1 - 

0 

1—1 

Nr  = 

Wr2  0 

,Ns  = 

Ws2  0 

0  I 

0  I 

(3.42) 


where 


W, 


Rl 


and 


now  denote  the  bases  of  the  null  spaces  of 


and  {Cy,Dye,Dyd),  respectively.  Provided  Wr2  and  Ws2  have  full  column  rank, 
this  leads  to  the  following  Generalized  Eigenvalue  Problem  (GEVP)  formulation 
[SLBB96] 

minimize  7  (3.43) 


subject  to 


R  I 
I  S 


>  0 


L3  G  L^,  J3  6  L^, 


L3  I 
I  J3 


>  0 


ZR  +  'y 

’  WIJ3WR2 

0 

’  WIJ3WR2 

0 

<  7 

0 

Js 

0 

J3 

Zs  +  7 

WJ,L3Ws2 

0 

'  IVJ2T3W52 

0 

<  7 

0 

L3 

0 

i3_ 

(3.44) 


(3.45) 


(3.46) 


(3.47) 


where  i?,  5,  T3,  J3,  J3,  £3,  .Bd,  Ge,  are  all  as  defined  in  Theorem  3.4.1. 

While  not  immediately  apparent,  the  matrix  sums  on  the  left-hand  sides  of  (3.46) 
and  (3.47),  when  carried  out,  become  independent  of  7.  Again,  this  problem  can 
be  solved  using  LMI  solvers  such  as  the  LMI  Control  Toolbox  in  MATLAB.  Often, 
Wr2  and  Ws2  do  not  have  full  column  rank;  when  this  occurs,  dummy  variables 
are  used  to  form  a  slightly  modified  GEVP  [GNLC92].  This  modified  GEVP  was 
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used  to  solve  for  the  F-18  controller  in  Chapter  IV.  Note  also  that,  depending  on 
the  problem,  the  number  of  dummy  variables  added  can  be  extremely  large.  This 
significantly  increases  the  computational  effort  required  to  solve  the  problem.  In 
general,  GEVP’s  require  more  computational  effort  than  testing  the  LMI  feasibility 
conditions,  and  adding  dummy  variables  further  increases  the  computer  time  and 
memory  requirements. 

Once  values  of  il,  5,  T3,  J3  have  been  found  for  some  feasible  performance  7 
(as  explained  in  Chapter  II,  the  chosen  solution  corresponds  to  the  analytic  center 
of  the  LMI),  the  LTI  controller  K  must  be  determined.  A  systematic  procedure  to 
derive  the  state-space  of  K  using  a  set  of  LMI  constraints  is  given  in  [AG95].  This 
method  offers  the  most  flexibility  since  additional  constraints  on  K  can  easily  be 
incorporated  in  the  form  of  more  LMI’s,  and  regularity  of  the  plant  is  not  required. 
Another  method  provides  a  series  of  analytical  expressions  to  calculate  the  individual 
state-space  matrices  of  K  [Gah94].  A  final  method  involves  simply  using  the  Riccati- 
based  central  Hoo  controller  equations  [DGKF89].  Regularity  of  the  plant  is  required 
for  this  last  approach;  however,  this  can  be  addressed  by  incorporating  small  weights 
on  all  the  measurements  and  control  inputs  of  the  system  to  ensure  regularity.  Since 
the  aircraft  design  system  developed  in  Chapter  IV  was  regular,  this  last  approach 
was  used  to  determine  the  controllers  in  this  thesis. 

3.5  Solving  the  Gain-Scheduled  Hco  Problem  for  Uncertain  LPV  Systems 

While  the  physical  parameters  9  are  assumed  to  be  measured,  there  may  be 
an  uncertainty  associated  with  the  measurements,  or  there  may  be  other  parameters 
which  remain  uncertain.  This  time-invariant  unmeasured  uncertainty  can  be  handled 
by  expanding  the  previous  results  to  make  use  of  classical  robust  control  techniques 
such  as  ^-synthesis. 

For  such  a  problem,  the  basic  LPV  structure  is  shown  in  Figure  3.4,  where  0* 
represents  the  block  operator  of  time-varing  parameters  as  defined  in  (3.5)  and  0„ 
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represents  the  constant  uncertainty  with  the  structure 

0„  =  blockdiag(^i/ri ,  •  •  • ,  ^5^5,  Ai, . . . ,  Af)  (3.48) 

A„  denotes  the  corresponding  structure  of  the  uncertainty  block 

A„  :=  {0„  ;  8i  €  C,  A,-  G  (3.49) 

0 

d 


Figure  3.4  LPV  control  structure  with  added  uncertainty 
The  total  “uncertainty”  structure  then  assumes  the  form 

Qt  0  0 

0  0t  0  (3.50) 

0  0  0„ 
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and  the  corresponding  set  of  scalings  is  then 


✓  r  “1  '' 

^  >  0  :  Lf  E  Lu  G  >  (3.51) 

0 

Figure  3.5  illustrates  the  transformed  LPV  control  structure,  taking  into  ac¬ 
count  the  added  time-invariant  uncertainty  block. 


Figure  3.5  Modified  LPV  structure  with  uncertainty 


Note  that,  in  order  to  perform  //-synthesis,  the  performance  channels  must  also 
be  wrapped  around  to  add  a  fictitious  0p  block  to  the  actual  uncertainty  block,  so 
that  all  of  the  various  objectives  can  be  considered  in  order  to  obtain  the  optimum 
//  scales.  This  is  shown  in  Figure  3.6.  The  corresponding  set  of  scalings  can  then  be 
written  as 


(3.52) 


0 

0 

_ 1 

0 

Lu 

0 

>  0  :  G  G  > 

1 - 

0 

0 

1 

a 

Figure  3.6  LPV / fi  control  structure 

The  D  scalings  defined  in  (3.52)  could  now  be  solved  for  in  place  of  the  L 
scalings  of  (3.26)  to  determine  the  7-suboptimal  gain- scheduled  controller.  Unfortu¬ 
nately,  this  augmented  scaling  resulting  from  the  addition  of  the  uncertainty  block 
0u  and  the  performance  block  0p  ruins  the  convexity  of  the  problem.  However, 
the  problem  is  convex  as  soon  as  the  Lu  and  Lp  blocks  of  the  scaling  matrix  are 
fixed,  since  we  then  essentially  restore  the  problem  defined  with  no  uncertainty.  In 
addition,  the  performance  7  can  be  optimized  using  ^-synthesis  when  Lt  is  fixed,  by 
performing  the  usual  D-K  iterations  [BDG'*'91].  This  motivates  a  mixed  LPV//X  ap- 


3-19 


proach  to  attempt  to  minimize  fi  and  optimize  7,  by  alternating  between  solving  the 
GEVP  and  performing  D-K  iterations.  Specifically,  the  following  ^^D-K-D”  scheme 
is  proposed  [SLBB96]: 


1.  Set  the  D  scale  defined  in  (3.52)  to  I. 

2.  Calculate  the  scaled  plant  Pa{s)  given  by 


D-^ 

0 


0  1 

> 

^r+w.2)  _  ^ 


(3.53) 


3.  Solve  the  LPV  problem  GEVP  for  Lt  using  the  scaled  plant  Pa- 

4.  Update  the  D  scale  with  Lt,  and  recalculate  Pa- 

5.  Design  an  LTI  controller  K  for  Pa{s)  using  either  the  Riccati-based  method  (if 
the  plant  is  regular)  or  the  LMI-based  methods  outlined  in  Section  3.4. 

6.  Fix  K  and  Lt,  and  find  scalars  me,  and  matrices  Mq^,Mq^  which  minimize 

rriQfl  0  0  0  0 

Ji=  0  M©„  0  Terf  0  M0„  0  ,  (3.54) 

0  0  Mq  0  0  Mq 

L  ^  J  L  J  00 

where  T^d  =  DFi{Pa,  K)D~^ .  This  is  identical  to  performing  the  first  D-K 
iteration.  Performing  additional  iterations  to  further  decrease  p,  is  not  rec¬ 
ommended,  since  this  can  quickly  result  in  high  order  Lu,Lp  scales,  thereby 
requiring  an  inordinately  large  computational  effort  to  solve  the  GEVP  LMI’s 
in  follow-on  iterations. 

7.  Normalize  m©,,  M©„,  and  M©p  by  m©,  and  obtain  and  Lp  by  fitting  transfer 
functions  to  the  normalized  M©’s.  Note  that,  for  m©,  to  be  scalar,  the  two  0^ 
parameter  blocks  must  be  combined  and  redefined  as  one  full  block  prior  to 
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step  6  [BDG+Ol].  This  is  necessary  in  order  to  be  able  to  normalize  the  Me  ma¬ 
trices,  although  it  does  unfortunately  introduce  substantial  conservativeness. 
Note  also  that  Lt  is  now  equal  to  1. 

8.  Fix  Lu  and  Lp  and  update  the  D  scale.  Restore  the  actual  uncertainty  block 
structure  (with  the  two  original  0t  scalar  blocks). 

9.  Repeat  steps  2  to  8  until  p,  no  longer  decreases. 

This  is  similar  to  //-synthesis  and,  likewise,  is  not  guaranteed  to  find  the  global 
minimum  solution;  however,  it  works  well  in  practice  [SLBB96],[BAG96]  and  will 
therefore  be  used  to  design  the  aircraft  controller  in  the  next  chapter.  Note  also 
that,  as  mentioned  previously,  an  alternative  to  solving  the  GEVP  is  to  iterate  on 
the  feasibility  problem  to  find  the  minimum  7  and  the  corresponding  T-scales.  This 
is  especially  recommended  when  the  number  of  dummy  variables  added  to  solve  the 
GEVP  is  extremely  large  and  would  thus  require  an  inordinate  amount  of  computer 
memory. 

Finally,  //-synthesis  is  strongly  recommended  even  when  uncertainty  is  not  con¬ 
sidered,  since  //  optimizes  the  scales  on  performance  and  robust  stability  to  further 
decrease  7;  whereas  the  LPV  approach  alone  will  optimize  only  the  scales  corre¬ 
sponding  to  the  time-varying  block.  In  fact,  the  standard  LPV  approach  is  simply 
steps  1  to  5  of  the  D-K-D  iteration  procedure. 
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IV.  Longitudinal  Aircraft  Control 

This  chapter  will  describe  the  implementation  of  gain  scheduling  theory  for 
LPV  systems  for  the  design  of  an  aircraft  pitch  rate  controller  operating  over  a  sub¬ 
stantial  portion  of  its  flight  envelope.  The  aircraft,  namely  an  F-18  Supermanuever- 
able  fighter,  will  be  designed  to  be  both  robustly  stable  and  meet  a  Level  1  perfor¬ 
mance  requirement  throughout  the  controller  design  envelope,  thereby  guaranteeing 
robust  performance. 

^.1  Longitudinal  Equations  of  Motion 

It  seems  worthwhile  to  first  briefly  review  the  derivations  and  assumptions 
involved  in  obtaining  the  aircraft  longitudinal  equations  of  motion  and  the  resulting 
state-space  dynamics. 

The  standard  nonlinear  equations  describing  motion  of  an  aircraft  are  given 
below.  Equation  (4.1)  describes  the  translational  forces  along  the  aircraft  body  axes, 
and  (4.2)  describes  the  rotational  forces  about  the  same  axes.  In  order  to  derive  these 
simplified  equations,  it  is  necessary  to  assume  that  the  aircraft  is  a  rigid  body,  that 
its  mass  remains  constant  over  time,  and  that  the  earth  provides  a  fixed  inertial 
reference  frame. 

X  =  m{U  +  QW  —  RV  +  g  sin  0) 

Y  =  m{V  -  PW  +  RU  -  gcosQsm^)  (4.1) 

Z  =  m(W  +  PV  —  QU  —  g  cos  0  cos  $) 

L  =  PIx  ~  RIxz  QR{lz  ~  ly)  ~  PQIxz 
M  =  QIy  +  PR{Ix  -  Q  -  R^Ixz  +  P^Ixz  (4.2) 

N  =  RIz  —  P Ixz  -|-  PQCy  —  Ix)  T  QRIxz 

In  the  above  equations,  U,  V,  W  are  the  translational  velocities  along  the  aircraft 
body  axes;  P,Q,R  are  the  rotational  velocities  about  the  same  axes;  Ix,Iy,Iz,Ixz 
are  the  moments  of  inertia;  g  is  gravity;  m  is  the  aircraft  mass;  and  X,  Y,  Z,  L,  M,  N 
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are  the  external  forces  and  moments  due  to  the  aircraft  aerodynamics  and  propulsion. 
The  Euler  angles  0,  W  describe  the  orientation  of  the  aircraft  with  respect  to  the 
earth 

0  =  <5  cos  $  —  i?  sin  $ 

$  =  P  +  Q  tan  0  sin  #  +  i?  tan  0  cos  $  (4-3) 

m  —  ^cos^  I  Qsin^ 
cos  ©  '  cos  0 

where  0  is  the  pitch  angle,  $  is  the  roll  angle,  and  is  the  yaw  angle. 

The  six  equations  of  motion  can  be  further  decoupled  into  two  sets  of  simulta¬ 
neous  equations:  three  equations  to  describe  the  longitudinal  dynamics,  and  three  to 
describe  the  lateral  dynamics.  By  assuming  that  the  aircraft  is  in  straight  and  level 
flight  with  only  longitudinal  disturbances  present,  we  can  neglect  all  lateral  forces 
and  moments,  resulting  in  the  following  simplified  equations: 

X  =  m{U  +  QW  +  g  sin  0) 

Z  =  m{W-QU  -gcosQ)  (4.4) 

M  =  QIy 

To  further  simplify  these  expressions,  it  is  then  assumed  that  the  velocities, 
forces,  moments  and  Euler  angles  can  be  considered  as  small  perturbations  from  their 
equilibrium  values;  and  that  products  of  these  perturbation  values  are  negligible. 
The  equations  are  then  expanded  to  express  the  longitudinal  forces  and  moments 
in  terms  of  changes  in  them  resulting  from  the  translational  and  angular  velocities 
perturbation  variables.  This  involves  determining  the  partial  derivatives  of  the  forces 
and  moments  with  respect  to  the  perturbation  variables.  A  detailed  development 
of  this  can  be  found  in  [Bla91]  or  any  good  text  on  aircraft  stability  and  control. 
The  resulting  linear  longitudinal  equations  of  motion  can  be  expressed  in  state  space 
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form  as: 
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where  state  variable  u  is  the  perturbation  velocity  along  the  x-axis,  a  is  the  pertur¬ 
bation  angle  of  attack,  q  is  the  perturbation  pitch  rate,  6  is  the  perturbation  pitch 
angle,  Uo  is  the  equilibrium  velocity,  and  0o  is  the  equilibrium  pitch  angle  (note: 
in  this  specific  context,  9  and  0o  indicate  pitch  angles,  not  time-varying  parame¬ 
ters).  The  A  matrix  is  composed  of  the  longitudinal  stability  derivatives,  while  the 
B  matrix  is  composed  of  the  longitudinal  control  derivatives.  For  convenience,  all 
these  derivatives  will  henceforth  be  referred  to  simply  as  the  stability  derivatives. 
The  control  input  is  8e  which  is  the  elevator  deflection.  Other  longitudinal  control 
inputs  are  possible  (such  as  trailing  edge  flaps  or  thrust  vectoring),  but  will  not  be 
considered  for  the  controller  design  in  this  thesis. 

Using  (4.1)-(4.2),  an  atmospheric  model,  and  a  precise  nonlinear  aircraft  model, 
local  equilibrium  flight  conditions  (also  called  trimmed  conditions)  can  be  determined 
where  all  accelerations  are  zero.  Since  this  could  allow  several  solutions  with  non¬ 
zero  velocities,  an  aircraft  flight  mode  must  also  be  specified;  the  most  common  of 
these  being  straight  and  level  flight.  Using  these  trimmed  flight  conditions,  values  for 
the  stability  derivatives  can  be  then  be  obtained.  In  this  way,  the  aircraft  dynamics 
at  a  point  in  its  flight  envelope  can  be  extremely  well  described  by  the  resulting  state 
space  equations. 

For  the  F-18  under  study,  such  trimmed  conditions  were  determined  at  18 
points  in  a  subsonic  flight  envelope  ranging  from  5000  ft  to  40000  ft,  with  Mach 
number  values  of  0.3  to  0.95,  using  the  nonlinear  F-18  aircraft  model  described  in 
[ABSB92].  These  points  are  shown  in  Figure  4.1,  with  the  corresponding  values  of 
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the  flight  conditions  given  in  Appendix  A  (the  state  space  A  and  B  matrices  are 
given  in  program  f  ISlongdat  .m  of  Appendix  C).  Note  that  JCq  and  Afg  are  often 
neglected  as  a  result  of  linearizing  the  equations  of  motion.  However,  these  stability 
derivatives  have  been  kept  in  this  case  since  the  trimmed  conditions  for  the  F-18 
provided  values  for  them. 


Mach  number 


Figure  4.1  Flight  envelope  showing  trimmed  data  points 

Typically,  an  aircraft  will  exhibit  two  natural  longitudinal  modes  of  motion: 
the  short  period  mode  and  the  phugoid  mode.  The  phugoid  mode  is  characterized 
by  a  long  time  period  (low  natural  frequency)  and  is  lightly  damped,  while  the 
short  period  mode  has  a  shorter  time  period  (higher  natural  frequency)  and  is  more 
heavily  damped.  For  a  manual  flight  control  system,  the  short  period  is  the  primary 
mode  of  interest  since  it  dominates  the  response  of  the  aircraft  to  pilot  inputs,  and 
the  pilot  can  usually  easily  correct  for  the  slight  error  caused  by  the  slowly  varying 
phugoid.  As  well,  the  pitch  rate  response  is  dominated  by  the  short  period  and,  for 
time  frames  of  typically  10  seconds  or  less,  exhibits  very  little  phugoid  motion.  For 
these  reasons,  the  longitudinal  state  space  model  can  be  reduced  to  a  second  order 
short  period  approximation,  given  in  (4.6).  Over  a  short  time  period,  it  provides  a 
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relatively  accurate  measure  of  the  aircraft  pitch  response. 
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In  this  thesis  two  pitch  rate  controllers  will  be  designed;  one  for  the  short 
period  model,  and  another  for  the  full  longitudinal  model.  However,  since  it  is  the 
more  accurate  one,  the  full  longitudinal  state  space  model  will  always  be  used  as  the 
aircraft  model  for  simulated  tests  of  the  controllers. 

4.2  The  LPV  Aircraft  Model 

As  explained  in  Chapter  III,  the  aircraft  must  first  be  modelled  as  an  LPV 
plant.  The  first  step  in  doing  this  is  to  model  each  of  the  stability  derivatives  as 
a  polynomial  expression  dependent  on  some  measurable  parameters  such  as  Mach 
number,  altitude,  angle  of  attack,  or  dynamic  pressure.  For  the  full  longitudinal 
model,  this  requires  curve-fitting  the  15  stability  derivatives  located  in  the  first 
three  rows  of  the  A  and  B  matrices  of  (4.5).  Note  that,  since  0o  is  unknown,  the 
expressions  in  the  last  column  of  the  A  matrix  are  also  being  treated  as  stability 
derivatives;  namely,  Xe  and  Zq  respectively. 

The  method  of  least-squares  was  used  to  fit  polynomial  expressions  of  various 
parameters.  Basically,  for  each  stability  derivative,  this  required  determining  which 
parameters  it  depended  upon  the  most,  and  combining  increasing  orders  of  these 
parameters  until  a  satisfactory  curve  fit  was  obtained.  While  not  all  of  the  stability 
derivatives  depended  as  well  on  some  parameters  as  others,  it  was  decided  that,  in 
this  case,  the  stability  derivatives  depended  mainly  on  Mach  number  and  altitude. 
These  were  therefore  chosen  as  the  time- varying  parameters  6{t)  on  which  the  LPV 
model  of  the  plant  will  depend.  Combining  various  powers  of  Mach  number  and 
altitude  to  obtain  as  good  a  fit  as  possible  while  maintaining  as  low  an  order  as 
possible  resulted  in  the  following  polynomial  functions  for  the  stability  derivatives. 
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Note  that  this  was  done  for  all  stability  derivatives,  although  for  convenience,  only 
the  short  period  stability  derivatives  are  listed  here. 


=  (8.4931e  +  00)  +  (-5.8946e-04)/i  +  (1.0695e-08)/i2  +  (-4.7974e  +  01)M 
+(3.0496e  -  0S)Mh  +  (-5.0945e  -  08)Mh^  +  (7.9105e  +  01)M^ 

+(-4.8436e  -  0S)M^h  +  (7.5817e  -  08)M^h'^  +  (-4.4981e  +  01)M^ 
+(2.5679e  -  03)M^h  +  (-3.7016e  -  08)M^h^ 

(4.7) 

Zg  =  (9.8313e  -  01)  +  (3.9110e  -  07)h  +  (-3.9162e  -  04)M  +  (-5.3330e  -  08)Mh 

(4.8) 

Zs,  =  (9.7204e  -  02)  +  (5.6449e  -  06)h  +  (-4.3877e  -  10)h^  +  (-1.0410e  +  00)M 
+(-5.5878e  -  06)Mh  +  (2.2222e  -  09)Mh^  +  (1.0255e  +  00)M^ 

+(2.2785e  -  05)M^h  +  (-4.1014e  -  09)M^k^  +  (-4.8758e  -  01)M^ 
+(-1.3154e  -  05)M^k  +  (2.3785e  -  09)MH^ 

(4.9) 

=  (4.0000e  +  02)  +  (-3.5864e  -  02)k  +  (8.1789e  -  07)k^  +  (-2.0713e  +  03)M 
+(1.7638e  -  01)Mh  +  (-3.9293e  -  06)Mk^  +  (3.4457e  +  03)M2 
+(-2.7712e  -  01)M^k  +  (6.0257e  -  06)M^h^  +  (-1.8966e  +  03)M^ 
+(1.4207e  -  01)MH  +  (-2.9944e  -  06)M^/i^ 

(4.10) 

Mg  =  (-1.9283e  +  00)  +  (1.0426e-04)h  +  (-1.8301e-09)k^  +  (7.3411e  +  00)M 
+(-2.3624e  -  04)Mk  +  (8.0762e  -  10)Mh^  +  (-1.4103e  +  01)M2 
+(1.7893e  -  04)M^  +  (8.9252e  -  09)M^h^  +  (6.7100e  +  00)M3 
+(3.4962e  -  05)M^h  +  (-9.5971e  -  09)M^h'^ 

(4.11) 

Ms,  =  (6.2187e  +  01)  +  (-5.6127e-03)/f  +  (1.0681e-07)/i2  +  (-3.i392e  +  02)M 
+(2.8099e  -  02)Mh  +  (-5.0915e  -  07)Mh^  +  (4.3654e  +  02)M2 
+(-4.1827e  -  02)M^h  +  (7.3379e  -  97)M'^h?  +  (-2.5310e  +  02)^^ 
+(2.1593e  -  92)M^h  +  (-3.4779e  -  97)M^h? 

(4.12) 
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Figures  comparing  the  polynomial  curve  fits  to  the  actual  values  for  all  15 
stability  derivatives  are  included  in  Appendix  B. 

As  explained  in  Chapter  III,  the  parameters  will  be  normalized  such  that  9{t)  < 

1  so  that  our  eventual  goal  will  be  to  design  a  controller  which  makes  the  robust 
performance  objective  7  <  1.  The  Mach  number  and  altitude  ranges  for  the  modelled 
envelope  are: 

M  G  [0.3,0.95]  and  h  G  [5000,40000]  ft.  (4.13) 

Expressing  these  in  terms  of  their  normalized  parameters  such  that  6{t)  G  [—1,1] 
gives: 

M  =  0.3250m  +  0.625  (4.14) 

h  =  175000,,  +  22500  (4.15) 

These  expressions  can  now  be  inserted  into  the  polynomial  functions  given  in  (4.7)- 
(4.12)  to  produce  normalized  polynomial  functions  of  the  stability  derivatives.  For 
Za,  this  would  result  in  the  following  equation: 

=  (-7.0152e  -  01)  +  (4.5365e  -  01)0/,  +  (-1.7379e  -  01)0^  +  (-4.8108e  -  O1)0m 
+(1.3941e  -  Ol)0M9h  +  (4.4512e  -  O1)0m0I  +  (-2.7859e  +  00)01^ 

+(4.8005e  +  00)01^0/,  +  (2.0740e  +  00)01^0^  +  (-2.0404e  +  00)0^ 

+(5.4195e  +  00)01^0/,  +  (-3.8915e  +  01)01^61 

(4.16) 

Once  all  the  equations  (4.7)-(4.12)  have  been  normalized,  we  now  have  a  non¬ 
linear  parameter- varying  model  of  the  plant,  where  each  of  the  elements  of  the  state 
space  model  is  given  by  a  polynomial  function  of  two  normalized  parameters:  0m 
and  Oh-  The  dependence  of  the  plant  on  these  parameters  is  known,  but  since  no 
information  on  Mach  number  and  altitude  are  known  a  priori  (except  their  range 
of  values),  this  fits  precisely  into  the  LPV  model  defined  in  (3.2).  We  now  need  to 
express  this  nonlinear  model  as  a  linear  fractional  function  of  the  parameters.  In 
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other  words,  we  need  to  model  this  as  a  nominal  plant  perturbed  by  the  parameters, 
with  added  elements  relating  how  the  parameters  will  affect  the  nominal  mapping. 

As  explained  in  Chapter  II,  an  LFT  can  be  thought  of  as  a  nominal  mapping 
perturbed  by  some  model  uncertainty,  with  the  other  elements  of  the  plant  relating 
how  the  uncertainty  affects  the  nominal  plant.  This  is  precisely  what  we  wish  to 
achieve,  although  in  this  case  the  “uncertainty”  is  composed  of  the  time-varying 
parameters.  Since  the  nominal  mapping  must  be  independent  of  the  parameters,  we 
can  clearly  see  that  the  nominal  plant  will  be  composed  of  the  constant  terms  in  the 
normalized  polynomial  expressions.  However,  in  order  to  determine  the  other  terms, 
it  is  first  necessary  to  convert  each  of  the  normalized  polynomial  equations  into  an 
LFT.  Looking  at  (4.5)  or  (4.6),  we  can  see  that  each  of  the  stability  derivatives  is 
simply  an  input/output  relation  between  specific  inputs  or  states  to  specific  outputs. 
Therefore,  we  can  use  the  technique  from  Chapter  II  to  convert  these  polynomial 
functions  to  LFT’s.  For  example,  the  short  period  equation  for  a  is 

a  =  ZaO!  -|-  ZgQ  Zs^Se  (4.17) 

where  Za  is  given  in  (4.16).  For  convenience,  we  can  rewrite  this  as 


d  =  Oioi  -|-  dg  -t-  d^j 


(4.18) 


where,  for  example, 


a. 


(4.19) 


The  equation  for  d^  is  now  a  simple  input/output  expression  whose  transfer  function 
is  given  by  the  Za  polynomial.  For  ease  of  notation,  (4.16)  will  be  generalized  as: 


Za  —  o,-\-h9h-\-c9\-\-d9  M-\-  f  ^ 

(4.20) 
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Following  the  procedure  outlined  in  Chapter  II,  the  block  diagram  shown  in 
Figure  4.2  can  now  be  drawn  to  determine  all  the  input/output  relations.  Note  that 
there  are  several  other  realizations  possible,  depending  on  how  the  0m’s  and  6hS 
are  arranged;  a  continuing  problem  in  control  theory  is  how  to  guarantee  that  we 
always  have  the  minimal  number  of  9  to  represent  polynomial  functions  such  as  these 
[ZDG96]. 


Figure  4.2  Block  diagram  of  Za  showing  input/output  relations 
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Pulling  out  the  0’s  into  one  block  diagonal  parameter  structure  and  writing 
the  outputs  in  terms  of  the  inputs  gives  (for  da): 
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or, 


Zc. 
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(4.22) 


where  is  a  coefficient  matrix  (not  a  transfer  matrix)  and,  as  expected,  the  constant 
term  =  a  is  the  nominal  map  for  the  upper  LFT  representation  of  Za-  For 
completeness. 


da  :=  ZaOL  =  F„(Za,0^<,)a,  where  0^^  = 


9m  h  0 

0  Ofih 


(4.23) 


This  LFT  representation  of  Za  is  shown  in  Figure  4.3.  Each  of  the  stability 
derivative  polynomials  can  be  converted  to  LFT’s  in  the  same  fashion.  When  this 
is  completed,  the  LFT’s  can  then  be  linearly  interconnected  to  form  the  parameter- 
varying  system  shown  in  Figure  4.4  (for  the  short  period  model). 
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This  linear  interconnection  of  LFT’s  can  now  be  grouped  into  one  combined 
LFT  by  gathering  the  known  systems  together  and  by  combining  the  varying  param¬ 
eter  blocks  together.  Using  the  notation  of  (4.23)  for  each  of  the  coefficient  matrices 
in  Figure  4.4,  we  can  obtain  the  following  state  space  equations  for  the  LPV  system. 
Note  that,  as  shown,  the  two  outputs  for  the  plant  have  been  chosen  to  be  the  states 
a  and  q. 
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(4.24) 


The  combined  parameter  block  0  has  a  structured  “uncertainty”  form  given  by 


0  =  hlockdiag{Qzcc^^Mai^Zg-,QMq,  ©z#,  ©mJ  (4.25) 
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Equivalently,  we  can  write  this  system  as  the  transfer  matrix  G(s)  perturbed 
by  the  parameter  structure  0.  G(s)  is  then 


D2 


0  C2 


0 


0  Dzs 


Dmcc  Dmo  0  0  Cm,  0  Cms  Dms 


BZa  0  0 


0  0  0  0 


5m.  0  0  Am.  0  0  0  0  0 


G{s)  = 


0  Bz,  0  0  A^,  0 


0  5, 


0  0  0  A, 
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0  0 


(4.26) 
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To  simplify  the  notation,  the  state  space  equations  will  be  grouped  together 
as  follows  _  _  _ 


Ag  Bg  Bxt,g 

=  Cg  Dgg  Dgug 

^yg  Bygg  DygUg 


(4.27) 


Therefore, 
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1 _ 

G{s)  =  Cg  Dgg  Dg^g  (4-28) 

Gyg  Bygg  DygUg 

Further  on,  we  will  want  to  extract  the  pitch  rate  output  q.  This  will  be  specified  as 

Q  =  Gyg^Xg.  (4.29) 


In  addition,  note  that  Dygg  =  0;  therefore,  assumption  3  of  Chapter  III  has  been 
met,  thereby  ensuring  that  the  LPV  controller  will  be  causal  and  well-posed. 


It  should  also  be  pointed  out  that,  while  the  ^m’s  and  O^s  for  each  individual 
stability  derivative  were  grouped  together  (as,  for  example,  in  (4.23))  these  repeated 
scalar  blocks  of  and  O^s  were  not  combined  together  in  (4.25)  for  the  different 
stability  derivatives.  For  practical  reasons,  it  is  usually  convenient  to  combine  them. 
This  then  requires  rearranging  the  corresponding  input/output  de,ee  channels  in 
(4.24)  and  (4.26).  This  will  be  assumed  to  have  been  performed  in  (4.27).  Thus  the 
time- varying  parameter  block  can  be  expressed  as 


^mItm  0 

0  dhlrh 


(4.30) 


where  vm  and  Vh  are  the  respective  dimensions  of  0m  and  0h-  Note  that  0  is  now  a 
repeated  scalar  structured  “uncertainty”  block. 

The  combined  LPV  aircraft  model  is  now  an  LFT,  and  can  be  simply  drawn 
as  in  Figure  4.5.  This  is  very  nearly  the  form  we  need  to  make  use  of  the  gain 
scheduling  theory  outlined  in  Chapter  III;  all  that  is  missing  is  the  overall  setup  of 
the  weighted  aircraft/controller  feedback  design  model. 


4.3  The  Design  Model 

The  problem  is  to  design  a  robustly  performing  pitch  rate  controller  for  the  F-18 
Supermaneuverable  fighter  aircraft.  Frequency  and  time  domain  specifications  for  a 
predicted  Level  1  handling  qualities  response  can  be  found  in  MIL-STD-1797A.  These 
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suggest  that  a  damping  factor  (sp  of  0.5  and  natural  frequency  Lo^p  of  4  radians/second 
would  meet  Level  1  handling  qualities  over  the  chosen  flight  envelope.  Thus  the  ideal 
pitch  rate  response  can  be  approximated  by  the  following  second  order  system: 


The  design  approach  will  be  to  model-match  the  response  of  the  plant  to  the 
reference  model  response  given  in  (4.31).  This  should  provide  nearly  the  same  pitch 
rate  response  over  the  entire  envelope.  Alternatively,  an  LPV  reference  model  could 
be  used  if  the  desired  response  had  to  vary  over  the  flight  envelope.  The  system 
design  model  is  shown  in  Figure  4.6. 


Figure  4.6  Design  model 
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In  state  space,  the  chosen  reference  model  can  be  written  as 


ih 
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.  . 

Converting  the  transfer  function  in  (4.31)  to  state  space  form  gives 
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(4.32) 


(4.33) 


Note  that  this  is  only  one  of  several  possible  state  space  forms  since  these  are  not 
unique. 

The  aircraft  plant  is  specified  in  (4.27)-(4.29).  Note  that,  as  shown  in  Figure 
4.6,  both  q  and  yg  are  output.  The  output  q  is  differenced  with  the  reference  pitch 
rate  response  qr  to  produce  an  error  signal  which  is  then  weighted  by  Wp  to  emphasize 
the  frequency  range  of  interest.  For  the  short  period  design,  Wp  was  chosen  as 

_  0.4(s  +  100) 


whereas  for  the  full  longitudinal  model,  Wp  was  chosen  to  be 


Wp  = 


0.04(s  +  100) 
5  +  4 


(4.35) 


For  the  short  period,  this  effectively  guarantees  that  the  steady-state  pitch  rate  error 
(to  a  unit  step  response)  in  the  final  design  will  be  less  than  0.1  degrees/ second  if 
an  Hoo  norm  (or  peak  y  value)  less  than  one  is  achieved.  It  can  also  be  thought  of 
as  a  low-pass  weight  on  sensitivity  to  provide  nominal  performance.  Some  iterating 
was  required  to  obtain  Wp,  because  the  aircraft  dynamics  and  the  design  problem 
formulation  determines  how  small  a  performance  error  will  be  possible  before  robust 
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performance  is  no  longer  met.  For  instance,  the  full  longitudinal  model  could  only 
guarantee  that  the  error  to  a  step  response  would  be  less  than  1  degree/second  for 
a  similar  flight  envelope.  This  is  only  a  conservative  guarantee  based  largely  on  the 
Small  Gain  Theorem  ;  in  fact,  the  response  will  generally  be  much  better.  The  state 
space  equations  for  Wp  can  be  generalized  as 

Xp  .Tp  -^p  Xp 

Cp  C'p  Dp  ^ 

The  output  yg  (composed  of  a  and  q)  is  added  to  a  noise  weighting  W„  to 
produce  a  measured  ym-  The  noise  weighting  can  be  thought  of  as  the  error  expected 
in  the  measurements  or,  alternatively,  as  a  weight  on  complementary  sensitivity.  This 
should  reduce  the  effect  of  the  high  frequency  measurement  noise  in  the  design.  In 
order  to  minimize  the  dimension  of  the  eventual  LMI  problem  and  the  number  of 
controller  states,  was  chosen  to  be  simply  a  static  weight  given  by 

Wn  =  0.05/2  (4.37) 

which  can  be  viewed  as  a  5%  error  expected  in  both  measurements,  or  as  a  weight 
on  complementary  sensitivity  between  the  output  yg  and  the  input  d„.  For  the  full 
longitudinal  design,  this  weight  was  decreased  to 

W„  =  0.01/2  (4.38) 

in  order  to  increase  the  relative  importance  of  performance  and  stability,  which,  for 
this  mode,  are  more  difficult  to  attain. 

Also  included  in  the  design  model  is  a  first  order  approximation  of  the  actuator 
dynamics  with  a  pole  at  s  =  —20.2.  As  shown  in  Figure  4.6,  we  would  like  not  only 
to  output  the  elevator  deflection  Sg,  but  also  the  rate  of  elevator  deflection  Sg-  This 


(4.36) 


allows  the  designer  to  penalize  either  the  elevator  usage  or  the  rate  of  elevator  usage 
(or  both).  The  state  space  equations  can  be  written  as 


Numerically,  this  leads  to 

Aq.  Ba 

Act{s)  = - 

D, 

Note  that,  in  this  case,  since  Se  =  Xa  then 

K  =  Xa=  AaXa  +  BaU  (4-41) 


=  (4.39) 

Ca  Da  U 


For  this  design,  penalizing  the  elevator  deflection  rate  was  sufficient.  The 
control  weight  or  penalty  Wc  was  therefore  selected  to  be 

Wc  =  l/70.  (4.42) 

This  ensures  that  the  maximum  allowable  deflection  rate  of  70  degrees/ second  will 
not  be  violated  if  an  Hoo  norm  (or  peak  jx  value)  less  than  one  is  achieved.  The 
control  output  Cg  can  then  be  written  as 

e,  =  WX.  (4.43) 


To  reflect  some  uncertainty  in  the  measured  parameters  (and  thus  yg)  and 
the  fact  that  we  are  not  really  modelling  high  frequency  dynamics,  an  uncertainty 
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weighting  Wu  was  added.  An  additive  uncertainty  “across”  G(s)  could  have  been 
used,  or  even  an  uncertainty  in  each  of  the  measured  time- varying  parameters,  but 
this  would  have  added  several  more  input /output  channels  and  increased  the  dimen¬ 
sion  of  the  LMI  problem  to  be  optimized.  Since  this  gain  scheduling  approach  is 
based  on  the  Small  Gain  Theorem  and  thus  quite  conservative  already,  it  was  finally 
decided  to  add  only  a  single  input  multiplicative  uncertainty.  After  some  iterations, 
the  weighting  was  chosen  as 

0.1(^  +  100) 

"  (s  +  10000) 

This  can  be  viewed  as  a  high-pass  weight  on  complementary  sensitivity  to  provide  ro¬ 
bust  stability,  or  as  increased  uncertainty  at  higher  frequencies.  It  can  be  generalized 
as 

Xu 

The  uncertainty  input  du  is  added  to  the  elevator  deflection  6e  to  produce  the  input 
Ug  of  the  aircraft  plant  G{s).  Specifically, 

Ug  =  Se  +  du.  (4.46) 

Finally,  to  increase  the  flexibility  of  the  controller,  a  2  degree-of-freedom  con¬ 
troller  design  was  adopted.  This  is  similar  to  the  stability  augmentation  system 
(SAS)  of  classical  control,  and  provides  the  controller  with  more  information  on  the 
states  a  and  q  than  it  would  have  with  only  a  standard  error  signal.  The  output  of 
the  controller  is  the  commanded  elevator  deflection  u  and  the  inputs  are  qc  and  ym, 
where  is  a  “noisy”  a  and  q  given  by 

ym=yg  +  Wndn.  (4.47) 


1 

e 

1 - 

1 _ 

Cu  Du 

- 1 

(0 

_ 1 

(4.45) 
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Therefore,  the  controller  input  y  is  defined  as 


y  =  (4.48) 

Vm 


The  linear  interconnection  of  all  these  components  as  depicted  in  Figure  4.6 
yields  the  following  state  space  equations  for  the  design  model: 
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We  can  now  write  this  as  the  design  plant  P{s)  given  by 


Pis)  = 
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(4.50) 


The  state  space  equations  can  then  be  grouped  together  to  recover  the  notation 
of  Chapter  III  as  follows 


A  Be  Bd  Bu 
Ce  Dee  Ded  Dou 
Ce  Dee  Ded  D  eu 
^yd  ^yu 


(4.61) 


where  e  =  gp  and  d  =  dp 


(4.52) 
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Therefore, 


’  A 

Be  Bd  Bu 

P(s)  = 

Ce 

Doe  Dod  Dqu 

Ce 

Dqo  D^d  Dqh^ 

_Cy 

^y6  Dyd  Dyy^ 

(4.53) 


This  is  now  exactly  the  LTI  design  plant  P{s)  referred  to  in  Chapter  III.  We 
can  now  add  the  time-varying  parameter  block  of  the  controller  and  the  resulting 
added  plant  interconnections  to  yield  the  augmented  plant  Pa(-s)  given  in  (3.21)  and 
reproduced  below. 
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(4.54) 


The  resulting  LPV  design  is  then  identical  to  the  one  shown  in  Figure  3.3  with  the 
time- varying  parameter  block  now  given  by 


A©  A  = 


0  0 
0  0 


(4.55) 


This  plant  Pa  is  the  one  used  to  solve  the  gain- scheduled  robust  performance 
Hoo  problem  described  in  Chapter  III.  By  solving  either  the  GEVP  or  iterating 
on  the  feasibility  problem,  one  can  determine  the  L-scales  required  to  find  the  cor¬ 
responding  controller.  Since  this  particular  design  includes  some  uncertainty,  the 
D-K-D  iteration  scheme  will  be  used  to  incorporate  the  advantages  of  //-synthesis. 
In  fact,  //-synthesis  is  recommended  even  when  uncertainty  is  not  considered,  since  // 
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optimizes  the  scales  on  all  input/output  channels;  whereas  the  LPV  approach  alone 
will  optimize  only  the  scales  corresponding  to  the  time-varying  block. 

Further,  for  the  purpose  of  this  thesis,  a  robust  stability  subproblem  will  be 
defined  where  the  only  input /output  channels  in  the  design  plant  are  dg  and  eg  and, 
correspondingly,  the  only  component  of  the  uncertainty  block  is  the  time-varying 
parameter  block  0.  This  can  also  be  achieved  by  making  all  the  design  weights 
extremely  small.  By  solving  the  GEVP  or  simply  testing  the  feasibility  problem  for 
7  <  1,  we  can  then  determine  if  the  plant  is  at  least  robustly  stabilizable  in  the  given 
flight  envelope.  If  not,  the  envelope  can  be  gradually  reduced  until  it  is.  Clearly, 
the  robust  performance  problem  will  require  an  envelope  at  least  as  small,  typically 
smaller. 

4.4  Implementation 

Most  of  the  programs  required  to  optimize  this  uncertain  gain-scheduled  Hoo 
problem  have  been  included  in  Appendix  C.  This  section  briefly  describes  these, 
along  with  some  other  implementation  issues. 

As  mentioned  previously,  the  full  longitudinal  state  space  A  and  B  matrices 
for  the  18  trimmed  data  points  are  available  in  program  f  ISlongdat  .m.  Using  these 
data  points,  the  program  flSpoly.m  performs  a  least-squares  curve  fit  to  provide  a 
polynomial  expression  for  each  of  the  15  stability  derivatives  (for  the  short  period, 
this  is  performed  by  f  ISsppoly  .m).  This  program  can  also  be  modified  to  change 
the  flight  envelope,  and/or  the  order  of  the  curve  fit,  and/or  the  desired  time- varying 
parameters.  For  the  simulation  aircraft  model,  the  full  longitudinal  model  was  used, 
and  its  stability  derivatives  were  fitted  with  polynomials  of  relatively  high  order 
(with  as  high  as  terms)  to  minimize  the  model  error.  Six  of  these  polynomials 

are  given  in  (4.7)-(4.12),  and  the  actual  curve  fits  are  illustrated  in  Appendix  B. 

Unfortunately,  these  high-order  fits  were  extremely  computationally  demand¬ 
ing  due  to  the  resulting  large  number  of  variables  to  optimize  when  the  LMI  problem 
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was  formulated,  especially  when  performing  D-K-D  iterations  which  add  additional 
weights  to  an  already  large  problem.  Computation  was  extremely  slow  and  often 
ran  out  of  available  memory.  In  order  to  improve  the  LMI  solver  efficiency,  the  or¬ 
der  of  the  polynomials  was  decreased  so  that  the  highest  terms  were  of  order  M'^h. 
In  addition,  after  several  design  iterations,  it  soon  became  clear  that  the  envelope 
was  too  large  for  only  one  controller.  For  this  reason,  the  envelope  was  iteratively 
decreased  until  robust  stability  was  at  least  attainable  (this  will  be  demonstrated 
in  Chapter  V),  and  then  further  decreased  until  robust  performance  was  possible. 
Thus,  the  design  envelope  for  both  the  short  period  design  and  the  full  longitudinal 
design  was  finally  chosen  to  be 

M  e  [0.5, 0.95]  and  h  G  [5000, 30000]  ft.  (4.56) 

The  resulting  curve  fits  for  these  envelopes  are  also  illustrated  in  Appendix 
B.  Due  to  the  reduced  envelope  (and  thus  less  points  to  curve  fit)  the  lower-order 
polynomials  do  model  the  aircraft  dynamics  relatively  well.  In  addition,  the  conser¬ 
vativeness  built  into  this  gain-scheduling  process  and  the  modelled  uncertainty  may 
in  large  part  offset  the  remaining  curve-fit  inaccuracies.  The  results,  as  we  will  see, 
tend  to  support  this.  Note  that  the  simulation  model  retains  the  entire  flight  enve¬ 
lope  and  the  original  higher  order  curve  fits  in  order  to  accurately  test  the  resulting 
controllers  beyond  their  chosen  design  envelopes. 

The  approach  of  the  last  section  was  then  used  to  normalize  the  polynomials 
and  convert  both  the  short  period  and  full  longitudinal  models  into  LFT’s.  As  you 
may  anticipate,  this  procedure  can  get  very  complex  and  intensive,  especially  consid¬ 
ering  the  fact  that  the  full  longitudinal  model  has  11  additional  stability  derivatives. 
Fortunately,  the  Wright  Laboratory  Flight  Dynamics  Group  has  developed  an  algo¬ 
rithm  to  automate  this  process  for  the  short  period  model  of  an  aircraft,  given  two- 
parameter  polynomial  expressions  for  the  stability  derivatives  [SLBB96].  A  slightly 
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modified  version  of  this  program  called  f  18sppoly21ft  .m  was  used  to  convert  the 
short  period  model  into  an  LFT.  This  program  was  then  adapted  to  convert  the  full 
longitudinal  aircraft  model  into  an  LFT;  it  is  included  as  program  f  18poly21ft  .m. 
In  addition,  the  programs  also  attempt  to  find  the  minimal  realization  of  the  result¬ 
ing  LFT.  As  mentioned  previously,  it  is  not  guaranteed  to  find  the  absolute  minimal 
realization,  but  it  does  manage  to  substantially  reduce  the  dimension  of  the  time- 
varying  parameter  block.  This  in  turn  reduces  the  number  of  input /output  channels 
and  thus  the  size  of  the  plant;  and,  as  a  result,  the  number  of  variables  which  the 
LMI  solver  will  eventually  be  attempting  to  optimize.  This  is  very  significant  since, 
as  mentioned,  the  LMI  solver  is  very  computationally  expensive  in  terms  of  both 
memory  and  CPU  requirements. 

Program  fl8ol.m  was  then  used  to  obtain  <5(5),  P{s),  and  Pa{s)  for  the  full 
longitudinal  model.  For  the  short  period,  this  is  performed  by  fl8olsp.m.  For 
this  design,  the  //-Toolbox  MATLAB  program  sysic.m  was  used  to  form  the  design 
plants;  this  automates  the  setup  of  the  plants  as  developed  in  the  last  section. 

We  now  have  the  design  model  necessary  to  perform  the  D-K-D  iteration, 
implemented  as  program  dkdM18  .mfor  both  design  models.  This  program  is  basically 
an  adaptation  of  the  //-Toolbox  dkit .  m  program  which  performs  the  D-K  iteration. 
Basically,  once  Po(-s)  has  been  created,  the  LMI  Control  Toolbox  is  used  to  solve 
either  the  GEVP  or  the  feasibility  problem  for  the  L-scales  corresponding  to  the 
time-varying  parameter  blocks.  For  a  large  problem,  such  as  the  full  longitudinal 
one,  it  is  often  faster  to  first  test  for  the  feasibility  of  desired  7  values;  this  helps 
to  quickly  eliminate  problem  formulations  that  will  not  lead  to  valid  designs.  The 
feasibility  problem  is  set  up  in  program  nshinf  lmi2  .m  to  test  the  LMI  constraints 
given  in  (3.36)-(3.39)  of  Chapter  III.  To  optimize  7,  or  to  solve  a  smaller  problem 
such  as  the  short  period  design,  the  GEVP  method  should  be  used.  Equations 
(3.36)-(3.47)  of  Chapter  III  are  used  to  set  up  the  GEVP  in  program  nshinf Imi.m 
and  solve  for  the  L-scales. 
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Specifically,  the  resulting  scales  are  the  Lfscales  of  step  3  in  the  D-K-D  itera¬ 
tion  procedure.  Since  Pa{s)  for  this  design  is  regular,  the  LTI  controller  can  then  be 
found  in  step  4  using  the  Riccati-based  Hoc  equations.  Recall  that  the  “uncertainty” 
block  at  this  point  is  given  by 

0i  0  0 

0  0t  0  (4.57) 

0  0  0„ 

where  each  0t  is  the  repeated  scalar  parameter  block  given  in  (4.30).  Since  we  will  be 
performing  //-synthesis,  we  must  add  a  fictitious  “uncertainty”  block  corresponding 
to  the  performance  channels.  Thus, 

0t  0  0  0 

0  0t  0  0  ^  ^ 

0LMI  =  (4.58) 

0  0  0„  0 

0  0  0  0p 

Prior  to  performing  the  //-synthesis  of  step  6,  the  two  0t  blocks  must  be  combined 
into  one  full  block.  Using  the  notation  of  (4.30),  Qt  would  then  become  a  square  full 
block  of  dimension  r2(M+h)-  Thus, 

0t  0  0 

0;.=  0  0„  0  (4.59) 

0  0  0p 

As  explained  in  Chapter  III,  this  is  required  in  order  for  to  be  scalar,  since 
m©,  must  be  used  to  normalize  the  M©  scales  corresponding  to  the  uncertainty  and 
performance  channels.  Performing  step  3  now  gives  the  scales  m©,,M©„,  and  Mep 
respectively.  Once  normalized,  the  scales  axe  fitted  to  transfer  functions  as  in  dkit  .m. 
The  resulting  Lu  and  Lp  scales  then  form  a  new  D-scale.  As  indicated  in  step  2,  Pa 
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is  then  scaled  to  form  Pa  and,  using  the  original  block  structure,  this  scaled  plant  is 
used  in  the  next  iteration  to  solve  the  GEVP  for  the  new  Lf-scales. 

This  process  is  repeated  in  dkdMlS.m  until  the  scales  converge  and  n  ceases  to 
decrease  or  is  less  than  one.  As  was  mentioned  in  Chapter  II,  /Lf-synthesis  is  not  yet 
implemented  for  real  or  even  scalar  blocks.  However,  we  can  perform  fi- analysis  on 
such  structures.  Since  our  time-varying  parameters  are  known  to  be  real,  we  can 
therefore  take  advantage  of  this  to  slightly  reduce  some  of  the  conservativeness  in  the 
process.  By  combining  and  reorganizing  the  two  original  Qt  blocks,  we  can  produce 
two  new  repeated  scalar  real  blocks  defined  as  follows 

0M  =  and  Qh  =  [Ohlrhl-  (4-60) 


Therefore, 


^Manalysis 


0Af  0 
0  Qk 
0  0 
0  0 


0  0 

0  0 

Qu  0 

0  0p 


(4.61) 


This  is  necessary  in  order  to  guarantee  that  each  identical  real  scalar  block  is  treated 
identically  (in  other  words,  by  grouping  all  the  Om's  and  6hS  together,  we  can  ensure 
that  they  all  represent  the  same  real  scalar  values  of  6m  and  9h).  Note  that  this  also 
requires  reorganizing  the  input/output  channels  of  the  scaled  plant  Pa.  While  no 
D-scales  are  generated,  this  does  produce  a  more  realistic  frequency  response  for  y. 


Even  though  the  actual  scales  used  are  found  using  0^  and  the  original  Pa,  as  the 
iterations  proceed  we  can  observe  the  more  realistic  value  of  y  decrease  as  well.  The 


peak  value  of  this  more  realistic  y  is  always  less  than  or  equal  to  that  of  the  complex 


y,  as  a.  result,  we  can  focus  on  this  more  realistic  response  and  stop  iterating  when 
it  decreases  below  one  (the  complex  y  at  this  point  is  usually  still  above  one). 
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While  there  is  no  guarantee  of  finding  the  optimum,  in  practice  the  program 
usually  converges  to  at  least  a  local  minimum  and,  as  will  be  demonstrated  in  Chap¬ 
ter  V,  there  is  a  substantial  drop  in  7  and,  therefore,  highly  improved  robust  perfo- 
mance  levels  compared  to  what  could  be  achieved  using  standard  Hoc  methods. 

On  occasion,  the  optimization  does  appear  to  diverge  for  a  few  iterations  before 
resuming  convergence.  This  may  be  due  to  the  presence  of  a  minimum  within  the 
iteration  tolerance.  One  of  the  LMI  solver  options  is  the  feasibility  radius,  and  if 
this  is  too  large,  the  resulting  scales  tend  to  be  discontinuous  and  impossible  to  fit 
transfer  functions  to.  This  also  leads  to  diverging  behavior.  If  the  radius  is  made 
too  small,  the  LMI  solver  may  not  find  an  existing  feasible  solution.  More  details  are 
given  in  [GNLC92].  As  well,  the  LMI  solver  occasionally  appears  to  search  endlessly 
for  a  feasible  solution  or  to  simply  “fall  asleep” ,  especially  when  trying  to  solve  very 
large  problems;  whether  this  is  due  to  the  LMI  solver,  MATLAB,  or  the  computer 
itself  is  unknown.  Thankfully,  this  unexplained  behavior  is  relatively  uncommon, 
and  only  manifests  itself  when  the  design  problem  grows  to  about  1500  variables  or 
more. 


4-28 


V.  Results  and  Simulations 


This  chapter  describes  the  approach  taken  to  define  the  design  envelope  for 
each  of  the  controllers,  and  describes  the  resulting  short  period  and  full  longitudinal 
controllers.  Finally,  several  simulations  will  be  performed  to  assess  their  actual 
performance  within  and  beyond  their  design  envelopes. 

5.1  The  Short  Period  Controller 

For  the  short  period  design,  testing  the  feasibility  of  the  robust  stability  sub¬ 
problem  indicated  that  a  stabilizing  LPV  controller  could  be  found  for  the  reduced 
flight  envelope  defined  by 

M  e  [0.5, 0.95]  and  h  e  [5000, 35000]  ft.  (5.1) 

Note  that  other  parts  of  the  flight  envelope  could  also  have  been  used  to  find  a 
robustly  stable  envelope;  the  decision  to  focus  on  this  portion  was  arbitrary.  For 
comparison,  the  Riccati-based  Hoo  method  was  used  to  check  the  feasibility  of  an 
LTI  controller  over  the  same  envelope;  the  resulting  norm  was  2.6938,  indicating 
that  no  LTI  controller  could  guarantee  stability. 

For  the  LPV  controller,  the  norm  of  the  gain-scheduled  Hoo  problem  was  found 
to  be  approximately  0.93.  While  stability  was  ensured,  this  left  very  little  room  in 
which  to  accomodate  the  performance  requirements  and,  as  a  result,  no  robustly 
performing  controller  could  be  found  which  met  the  objective  7  <  1.  Since  the 
problem  formulation  is  quite  conservative,  the  solution  can  be  expected  to  withstand 
larger  deviations  in  M  and  h  than  it  was  actually  designed  for.  For  this  reason,  rather 
than  decrease  the  weight  on  performance,  the  size  of  the  design  envelope  was  slightly 
reduced  to 

M  G  [0.5, 0.95]  and  h  6  [5000, 30000]  ft.  (5.2) 
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Using  this  envelope,  the  plant  Pa{^)  was  generated  as  shown  in  Chapter  IV. 
A  check  of  the  Hoo  norm  of  the  unsealed  problem  gave  7  =  2.7325,  which  exceeds 
the  desired  objective.  (Also,  testing  for  a  robustly  stabilizing  LTI  controller  for  this 
smaller  envelope  still  gives  an  excessive  norm  of  2.2365.)  Performing  the  D-K-D 
iterations  gives  the  data  shown  in  Table  5.1 


Table  5.1  Short  Period  Design  D-K-D  Data 


Iteration 

# 

Number  of 
variables 

Hoo 

norm 

Controller 

order 

Peak 

value 

P eak  ^analysis 
value 

1 

480 

1.5777 

7 

1.4920 

1.4075 

2 

626 

1.1301 

13 

1.1364 

1.0647 

3 

824 

1.0116^ 

19 

1.0212 

0.9638 

4 

1436 

0.9826 

31 

0.9805 

0.9136 

Note  that  the  solution  to  the  standard  gain-scheduled  Hoo  problem  (without 
doing  any  /;i-synthesis)  is  given  in  iteration  1  and  would  only  reach  7  =  1.5777,  which 
points  out  the  advantages  of  performing  the  D-K-D  procedure. 

The  controller  generated  for  either  of  iterations  3  or  4  can  be  used,  since  both 
indicate  that  the  true  value  of  //  is  less  than  1;  however,  the  controller  of  iteration 
3  is  selected  since  it  has  a  much  lower  order,  further  emphasizing  the  benefits  of 
performing  the  /^analysis  step.  Even  so,  the  order  of  the  controllers  can  be  quite 
high  depending  on  the  order  of  the  plant,  the  weights,  and  especially  the  order  of 
the  rational  D-scales  which  are  applied  to  each  of  the  channels.  For  this  design, 
the  plant  and  weights  contribute  7  states,  whereas  the  Lp  portion  of  the  D-scale 
contributes  12  states  (2nd  order  rational  functions  were  used  on  the  performance 
channels).  As  a  result,  the  controller  has  19  states.  This  is  relatively  high,  but 
reasonable  considering  this  single  controller  can  operate  over  a  large  portion  of  the 
envelope.  The  order  can  also  be  reduced  via  model-order  reduction,  but  that  was 
not  done  for  this  thesis.  Finally,  it  should  be  pointed  out  that  the  controller  order 
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does  not  necessarily  grow  with  each  iteration;  it  only  does  so  in  this  case  since  larger 
Z)-scales  were  required  at  each  iteration  to  better  approximate  and  decrement  /x. 

Due  to  the  relatively  high  order  of  the  controller,  frequency  response  plots  of 
the  controller  can  be  used  to  demonstrate  that  the  resulting  LPV  controller  does 
indeed  adapt  to  changing  values  of  M  and  h.  The  particular  plots  shown  in  Figure 
5.1  represent  the  frequency  response  of  the  controller  output  u  with  respect  to  the 
commanded  pitch  rate  qc  for  various  values  of  the  parameters. 


Figure  5.1  Frequency  response  of  the  short  period  LPV  controller 


5.2  The  Full  Longitudinal  Controller 

For  the  full  longitudinal  design,  the  feasibility  of  the  robust  stability  subprob¬ 
lem  was  first  tested  to  determine  the  envelope  over  which  a  stabilizing  controller 
could  be  found.  This  resulted  in  a  slightly  smaller  envelope  than  the  one  found  for 
the  short  period  design,  indicating  that  the  phugoid  mode  does  increase  the  potential 
for  instability.  The  stabilizable  envelope  chosen  is  defined  by 

M  e  [0.5, 0.95]  and  h  e  [5000, 30000]  ft.  (5.3) 
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For  comparison,  the  feasibility  of  an  LTI  controller  over  the  same  envelope  was  again 
tested;  the  resulting  Hoc  norm  was  157.41,  indicating  that  no  LTI  controller  could 
even  come  close  to  guaranteeing  robust  stability. 

For  the  LPV  controller,  the  norm  of  the  gain-scheduled  Hoo  problem  was  found 
to  be  approximately  0.84,  thereby  guaranteeing  robust  stability.  In  this  case,  de¬ 
creasing  the  envelope  to  meet  the  robust  performance  objective  was  not  sufficient; 
in  fact,  no  sizable  envelope  was  found  which  could  provide  robust  performance  using 
the  same  weights  as  the  short  period  design.  Also,  for  very  small  envelopes,  the 
LPV  plant  becomes  based  on  only  a  few  trimmed  points;  in  combination  with  the 
conservatism  of  the  method,  the  results  then  become  questionable.  A  better  ap¬ 
proach  at  that  point  may  simply  be  to  model  the  plant  as  LTI  and  add  some  general 
uncertainties,  as  in  a  standard  Hooln  design. 

In  any  event,  in  order  to  attempt  to  gain-schedule  the  full  longitudinal  model, 
the  noise  weight  Wn  and  the  performance  weight  Wp  were  decreased  to  the  levels 
described  in  Chapter  IV  in  order  to  direct  more  controller  effort  towards  maintain¬ 
ing  robust  stability.  Of  course,  as  was  mentioned,  in  doing  so  we  are  allowing  some 
degraded  performance  in  the  form  of  greater  error  and  more  susceptibility  to  noise. 
After  some  initial  trials,  it  was  decided  to  attempt  to  gain-schedule  the  entire  stabi- 
lizable  envelope,  since  slightly  smaller  envelopes  required  just  as  large  a  decrease  in 
the  performance  weight. 

Using  this  envelope,  the  plant  Pa{s)  was  then  generated  for  the  full  longitudinal 
design  problem.  A  check  of  the  Hoo  norm  of  the  unsealed  problem  gave  7  =  157.43, 
which  far  exceeds  the  desired  objective.  (Recall  that  the  robustly  stabilizing  LTI 
controller  for  this  envelope  gives  a  norm  of  157.41).  Performing  the  D-K-D  iterations 
gives  the  data  shown  in  Table  5.2 

As  before,  the  solution  to  the  standard  gain-scheduled  Hoo  problem  is  given  in 
iteration  1  and  reaches  7  =  1.1934.  The  controller  generated  for  any  of  the  last  three 
iterations  can  be  used,  since  they  all  indicate  that  the  true  value  of  //  is  less  than  1; 
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Table  5.2  Full  Longitudinal  Design  D-K-D  Data 


Iteration 

# 

Number  of 
vaxiables 

Hoo 

norm 

Controller 

order 

Peak  n 
value 

P eak  //analysis 

value 

1 

934 

1.1934 

9 

1.1934 

1.0377 

2 

1178 

1.0153 

17 

1.0199 

0.9976 

3 

1112 

1.0146 

15 

0.9966 

4 

1112 

1.0155 

15 

Dm 

0.9919 

however,  the  controller  of  iteration  4  is  selected  since  it  has  a  lower  order  and  slightly 
lower  fj,  value.  In  this  case,  the  plant  and  weights  contribute  9  states,  whereas  the 
Lp  portion  of  the  D-scale  contributes  3  states  (1st  order  rational  functions  were  used 
on  the  performance  channels).  As  a  result,  the  controller  has  15  states.  The  order 
can  then  be  further  reduced  via  model-order  reduction.  Note  that  the  order  of  this 
controller  happens  to  be  lower  than  that  for  the  short  period  controller.  This  is  due 
to  the  D-scales  which  are  generated  as  rational  approximations  to  the  actual  scales. 
These  in  turn  depend  on  how  “smooth”  the  LMI  and  //synthesis  scales  turn  out  to 
be  as  a  result  of  the  problem  formulation,  tolerances,  feasibility  radius  of  the  LMI 
solver,  and  numerical  conditioning  of  the  problem  to  be  solved. 

Once  again,  frequency  response  plots  of  the  controller  can  be  used  to  demon¬ 
strate  that  the  resulting  LPV  controller  does  indeed  adapt  to  changing  values  of  M 
and  h.  The  particular  plots  shown  in  Figure  5.2  represent  the  frequency  response  of 
the  controller  output  u  {SeJ  with  respect  to  the  commanded  pitch  rate  qc  for  various 
values  of  the  parameters. 

5.S  Simulations 

Several  simulations  were  carried  out  to  validate  both  the  short  period  and  full 
longitudinal  controller  designs.  Unfortunately,  the  nonlinear  aircraft/atmospheric 
model  was  not  available  to  perform  true  dynamic  simulations  of  the  controllers; 
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Figure  5.2  Frequency  response  of  the  full  longitudinal  LPV  controller 

therefore,  MATLAB’s  simulation  toolbox  SIMULINK  [SIM]  was  used.  Two  types 
of  simulations  were  carried  out;  the  first  type  is  “static”  and  tests  the  controllers  at 
arbitrary  points  within  and  outside  their  design  envelopes  using  LTI  aircraft  models 
accurate  at  the  selected  test  points.  For  small  changes  in  pitch  rate  (which  is  an 
assumption  of  the  linearized  equations  of  motion),  this  should  provide  a  realistic 
transient  response  as  long  as  the  state  and  input  variables  remain  relatively  small. 
The  second  type  of  simulation  attempted  is  “dynamic”,  and  tries  to  account  for  the 
changing  aircraft  model  as  the  aircraft  moves  through  the  flight  envelope.  While  this 
is  potentially  more  realistic,  it  is  still  nevertheless  limited  by  the  assumptions  that 
were  made  in  order  to  linearize  the  equations  of  motion.  Finally,  simulations  are 
also  carried  out  for  LTI  controllers  designed  using  standard  /[^-synthesis  methods,  in 
order  to  compare  their  behavior  to  that  of  the  LPV  controllers. 

5.3.1  The  “Static”  Controller  Simulations.  Figure  5.3  illustrates  the 
SIMULINK  setup  used  to  test  the  controllers  at  arbitrary  test  points  in  the  flight 
envelope.  At  such  test  points,  the  values  of  M  and  k  are  known  and  thus  the  LTI 
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aircraft  model  can  be  generated  by  performing  an  upper  LFT  transformation  of  the 
high-order  full  longitudinal  simulation  LPV  aircraft  model  with  the  known  param¬ 
eter  block.  A  similar  lower  LFT  transformation  is  performed  on  the  controllers  to 
convert  them  from  LPV  controllers  to  LTI  controllers  tuned  to  the  test  point.  As 
illustrated,  the  simulation  also  includes  wind  and  measurement  noises,  as  well  as 
rate  and  saturation  limiters  on  elevator  deflection. 

For  the  static  short  period  controller,  the  first  four  simulations  performed  cor¬ 
respond  to  the  corners  of  the  design  envelope  formed  by  the  four  combinations  of 
minimum  and  maximum  values  of  the  parameters  M  and  h.  For  clarity,  these  step 
responses  are  simulated  without  the  added  noises.  The  results  are  shown  in  Figure 
5.4.  A  fifth  test  uses  the  parameter  values  at  the  center  of  the  design  envelope  and 
adds  the  noises.  The  results  for  this  simulation  are  shown  in  Figure  5.5. 

The  pitch-rate  response  displays  predicted  Level  1  flying  qualities  with  excellent 
tracking,  minimal  overshoot  and  no  oscillatory  behavior.  Even  though  the  aircraft 
simulation  model  is  full  longitudinal,  the  phugoid  mode  does  not  noticeably  affect 
the  time  response.  These  first  five  tests  indicate  that  the  controller  should  work  well 
“statically”  in  at  least  the  design  envelope. 

Three  additional  simulations  were  performed  at  arbitrary  points  outside  the 
design  envelope  to  examine  the  behavior  of  the  controllers  over  areas  in  which  robust 
performance  is  no  longer  guaranteed.  For  clarity,  these  tests  were  performed  without 
the  added  noises.  The  results  are  shown  in  Figure  5.6.  While  all  three  cases  are  out¬ 
side  the  robust  performance  envelope,  note  that  the  first  one  (M  =  0.45,  h  =  10000 
ft)  is  still  within  the  robust  stability  envelope  and  does  quite  well,  an  indication  that 
the  robust  performance  envelope  is  conservative.  The  next  case  {M  =  0.35,  h  =  5000 
ft)  is  not  as  good,  but  still  indicates  that  the  robust  stability  envelope  is  conserva¬ 
tive.  Finally,  the  third  case  (M  =  0.3,  h  =  4000  ft)  is  too  far  beyond  the  stability 
envelope  for  the  controller  to  even  maintain  stability  (because  of  this,  only  a  few 
seconds  are  shown  on  the  plots). 
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Figure  5.5 


Static  short  period  controller  simulation  at  the  center  of  the  design 
envelope  (with  added  noises) 


5-10 


q  command  (deg/sec) 


time  (sec) 


time  (sec) 


Figure  5.6  Static  short  period  controller  simulations  beyond  the  design  envelope 
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In  the  case  of  the  full  longitudinal  controller  design,  the  results  are  far  dif¬ 
ferent  and  are  illustrated  in  Figure  5.7.  By  choosing  to  gain-schedule  a  fairly  large 
envelope,  we  have  had  to  decrease  the  weight  on  performance  in  order  to  main¬ 
tain  stability.  Thus,  the  results  are  disappointing  but  not  unexpected.  The  first 
three  cases  shown  are  all  at  the  edge  or  within  the  design  envelope,  and  stability 
is  indeed  maintained  at  the  expense  of  any  visible  performance.  As  well,  the  last 
case  {M  =  0.45,  h  =  10000  ft)  illustrates  how  easily  stability  is  lost  just  outside 
the  design  envelope.  As  a  result  of  these  simulations,  other  designs  were  attempted 
with  smaller  envelopes;  unfortunately,  until  the  envelope  becomes  extremely  small, 
the  results  are  not  much  better.  In  retrospect,  this  appears  to  be  largely  due  to 
the  9  additional  stability  derivatives  and  the  correspondingly  larger  parameter  block 
which  undoubtedly  increase  the  conservativeness  of  the  problem.  Specifically,  for 
the  previous  short  period  design,  the  varying  parameter  blocks  for  the  plant  and 
controller  each  consisted  of  3  ^m’s  and  3  0hS,  while  for  the  full  longitudinal  design 
they  consisted  of  5  ^m’s  and  5  6hS,  even  though  both  were  curve-fitted  with  poly¬ 
nomials  of  order  M^h.  Due  to  the  Small  Gain  Theorem,  these  parameters  are  then 
treated  as  complex  variables;  this  is  not  only  very  conservative,  but  due  to  the  the 
greater  number  of  them  in  the  full  longitudinal  problem,  this  effectively  adds  much 
more  demanding  constraints  to  the  problem.  The  LTI  norms  for  the  robust  stability 
subproblems  are  a  good  indication  of  this;  for  the  full  longitudinal  design  the  norm 
was  70  times  greater  than  the  one  for  the  short  period  design,  indicating  how  much 
more  difficult  the  full  longitudinal  design  is  to  control.  Related  to  this  is  the  fact 
that,  in  //-synthesis  alone,  the  parameter  block  forms  a  much  larger  complex  full 
block  than  was  the  case  for  the  short  period.  Unfortunately,  there  is  little  that  can 
be  done  to  avoid  this  in  the  current  method,  except  trying  to  keep  the  order  of  the 
curve-fits  to  a  minimum,  and  having  as  few  parameter-dependent  elements  of  the 
plant  as  possible. 
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Finally,  it  should  also  be  pointed  out  that  the  ideal  2nd  order  reference  model 
does  not  make  any  allowances  for  the  phugoid  mode.  A  more  complex  Level  1  per¬ 
formance  reference  model  (or  perhaps  a  parameter- varying  reference  model)  might 
allow  a  greater  weight  on  performance  and  thus  improve  the  resulting  level  of  per¬ 
formance. 

5.3.2  The  “Dynamic”  Controller  Simulations.  For  more  realistic  simu¬ 
lations,  both  the  aircraft  model  and  the  controller  should  remain  LPV  and  thus 
depend  upon  the  changing  operating  conditions.  As  in  a  real  implementation,  this 
then  requires  updating  the  controller  and  the  aircraft  with  the  latest  available  values 
of  the  M  and  h  parameters  so  that  the  LPV  controller  and  the  aircraft  can  adapt  to 
the  changing  flight  envelope.  However,  since  M  and  h  are  not  directly  available  from 
the  “outside  world”  or  a  precise  non-linear  model  of  the  aircraft/atmosphere,  we  can 
approximate  them  using  the  states  of  the  aircraft  model.  This  can  be  achieved  at 
each  simulation  time-step  i  using  the  following  relationships  [Bla91],[BS89]: 

Ui  =  Ui-i  -^  Ui 

Ah, 

h{  —  h{ — 1  -j-  ^^hi 

Ti  =  518.67  -  0.003565/ii 

Ui 

1.4  X  1716.16Ti 

where  U  is  the  total  aircraft  velocity,  u  is  the  perturbation  velocity,  0  is  the  per¬ 
turbation  pitch  angle,  a  is  the  perturbation  angle  of  attack.  Ah  is  the  change  in 
altitude,  and  T  is  the  ambient  temperature  in  degrees  Rankin.  Note  that  u,  0,  and  a 
are  available  as  state  variables.  Since  the  trimmed  conditions  were  given  in  degrees, 
a  conversion  to  radians  is  required  to  obtain  Ah.  Also,  by  choosing  some  Mq  and 
ho  as  the  initial  position  in  the  flight  envelope,  we  can  calculate  Uq  by  using  the  last 
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Figure  5.7  Static  full  longitudinal  controller  simulations 


two  equations.  Finally,  the  equation  for  T,  is  only  valid  for  up  to  36000  feet  altitude; 
this  is  sufficient  for  all  the  test  cases. 

In  a  linear  simulation,  the  state  space  models  of  the  controller  and  plant  are 
generally  fixed,  while  the  states  and  input/output  channels  vary  with  time.  This 
leads  naturally  to  think  of  M  and  h  as  the  output  of  a  block  performing  the  functions 
given  in  (5.4),  and  whose  inputs  are  the  states  u,  0,  and  a.  The  outputs  would  then 
be  normalized  and,  since  they  form  the  parameter  block,  be  fed  back  to  both  the 
aircraft  model  and  the  controller  model.  An  example  of  this  for  the  short  period 
controller  is  illustrated  in  Figures  5.8  and  5.9.  Unfortunately,  this  does  not  work  in 
SIMULINK  due  to  the  unusually  high  number  of  algebraic  loops  which  are  created 
and  would  have  to  be  resolved  at  each  time-step  iteration. 

Fortunately,  SIMULINK  allows  the  user  to  define  his  own  model  blocks  using 
what  are  called  S-functions.  This  technique  was  used  to  define  two  new  state-space 
blocks  whose  state-space  matricies  could  be  updated  at  each  iteration.  These  new 
blocks,  called  LPVG  and  LPVK,  are  used  to  update  the  aircraft  model  and  controller, 
respectively,  at  each  time-step  in  the  simulation.  The  details  of  the  new  state- 
space  blocks  are  given  in  programs  LPVG.m  and  LPVK.m  of  Appendix  D.  The  final 
SIMULINK  setup  is  given  in  Figure  5.10  (the  uat2Mh  functional  block  is  identical  to 
that  in  Figure  5.9). 

Two  tests  using  the  short  period  controller  will  be  performed  using  this  dy¬ 
namic  simulation.  It  should  be  emphasized  that  this  is  only  an  approximation  to 
a  full  nonlinear  aircraft/atmospheric  model  and,  as  such,  is  not  meant  to  exactly 
duplicate  the  behavior  of  the  aircraft.  There  are  several  aircraft  attributes  (such  as 
added  drag  and  loss  of  lift)  which  are  lost  when  only  steady  level  flight  conditions 
form  the  basis  of  a  simulation  model.  Still,  the  aircraft  simulation  model  appears 
to  behave  appropriately  when  compared  to  actual  published  flight  trajectories  of 
similar  fighter  aircraft. 
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Figure  5.8  Possible  dynamic  LPV  simulation  model 
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Figure  5.9  uat2Mh  functional  block  (converts  states  u,  a,  6  to  parameters  M,  h) 
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CO 


Figure  5.10  Dynamic  LPV  simulation  model 
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Several  pitch-rate  commands  were  chosen  to  test  as  much  of  the  envelope  as 
possible  without  having  access  to  a  thrust  input.  The  resulting  dynamic  simulations 
are  illustrated  in  Figures  5.11  to  5.22.  The  first  four  figures  illustrate  how  the 
controller  reacts  to  a  doublet  pitch-rate  command;  while  the  next  four  illustrate  how 
it  reacts  to  the  opposite  pitch-rate  command  with  added  noises.  The  last  four  figures 
test  the  uppermost  part  of  the  envelope  using  basic  pitch-rate  command  pulses.  The 
results  for  these  dynamic  tests  are  once  again  quite  good  with  excellent  tracking 
and  noise  response,  indicating  that  Level  1  flying  qualities  can  be  expected  using 
this  controller  within  the  design  envelope.  As  well,  even  though  the  commands 
are  held  for  relatively  long  periods  of  time,  the  phugoid  has  no  noticeable  effect  on 
the  response.  This  suggests  that  this  short  period  design  may  well  be  sufficient  for 
adequate  pitch-rate  control,  and  that  a  full  longitudinal  design  may  not  be  necessary. 

A  Comparison  with  LTI  fi-Synthesis  Controllers 

It  has  already  been  demonstrated  in  Chapter  IV  that  the  LPV  controller  sta¬ 
bilizes  a  much  larger  envelope  than  any  LTI  controller  could.  It  would  also  be 
interesting  to  compare  the  performance  of  optimal  LTI  controllers  to  that  of  the 
LPV  controllers.  The  best  way  to  do  this  is  by  examining  their  H^o  norms  and  time 
responses. 

Using  exactly  the  same  weights  and  LPV  design  plant  P{s)  developed  in  Chap¬ 
ter  IV  and  specific  parameter  values  of  M  and  h,  we  can  obtain  a  plant  Pf  by  the 
upper  LFT  transformation 

P;(s)  =  F„(P,0).  (5.5) 

This  is  now  an  LTI  plant,  and  standard  Riccati-based  methods  can  then  be 
used  to  find  an  optimal  Hoo  controller.  For  a  fair  comparison,  the  LTI  design  model 
should  also  include  some  model  uncertainty  in  order  to  account  for  the  varying  nature 
of  the  plant;  otherwise,  it  will  have  excellent  performance  only  at  the  design  point. 
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Figure  5.14  Dynamic  response  to  a  positive  doublet  pitch-rate  command:  a,  6 
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Dynamic  response  to  a  negative  doublet  pitch-rate  command  with 
added  noises:  qc,  q 
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Figure  5.17  Dynamic  response  to  a  negative  doublet  pitch- rate  command  with 
added  noises:  Se,u 
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Figure  5.18  Dynamic  response  to  a  negative  doublet  pitch-rate  command  with 
added  noises:  a,  0 


5-27 


q  (deg/sec)  q  command  (deg/sec) 


20 


10 

0 

-10 


(I _ I _ 1 _ ^ ^ - 1 - 1 

0  5  10  15  20  25  30 

time  (sec) 


5-28 


5-29 


^  ^  elevator  (deg) 


5-3i 


pitch  angle  (deg)  aoa  (deg) 


Figure  5.22  Dynamic  response  to  pitch-rate  command  pulses;  a,  9 
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The  size  of  the  uncertainty  chosen  will  affect  the  tradeoff  between  performance  and 
stability;  as  a  result,  several  iterative  designs  should  be  performed  until  the  right 
mix  is  achieved.  Designs  along  these  lines  revealed  that  no  LTI  controller  could 
guarantee  robust  performance  for  a  sizable  portion  of  the  design  envelope;  therefore, 
for  the  purpose  of  this  comparison,  the  multiplicative  uncertainty  already  in  the 
design  model  was  used  in  combination  with  ^-synthesis  to  yield  several  controllers  of 
different  orders  optimized  for  the  center  of  the  LPV  design  envelopes.  The  controllers 
which  gave  the  most  acceptable  performance  throughout  the  LPV  design  envelope 
when  simulated  against  the  full  longitudinal  simulation  model  were  then  used  for 
comparisons. 

For  the  short  period  design,  the  resulting  robust  performance  controller  re¬ 
sulted  in  a  peak  n  value  of  0.326,  while  the  nominal  performance  norm  was 
0.2076;  both  indicate  that  robust  performance  for  a  sizable  envelope  should  not  be 
difficult  to  achieve  and  that  a  larger  uncertainty  could  be  addressed.  On  the  other 
hand,  the  full  longitudinal  design  had  a  peak  value  of  0.982  and  nominal  perfor¬ 
mance  Hoo  norm  of  0.986;  this  emphasizes  the  fact  that,  using  the  desired  reference 
model,  robust  performance  will  likely  not  be  achievable,  since  there  is  very  little 
room  to  satisfy  both  nominal  performance  and  robust  stability  for  even  the  small 
multiplicative  uncertainty.  This  was  demonstrated  in  the  full  longitudinal  LPV  de¬ 
sign,  where  we  could  find  a  controller  to  maintain  stability  only  at  the  expense  of 
performance. 

Static  step  responses  for  the  LTI  controllers  when  tested  on  the  full  longitu¬ 
dinal  simulation  model  are  shown  in  Figures  5.23  and  5.24.  As  expected,  the  full 
longitudinal  LTI  controller,  which  is  based  on  a  plant  similar  to  that  of  the  simula¬ 
tion  model,  does  well  at  the  design  point  (M  =  0.75,  h  =  15000  ft)  but  off-design 
points  show  degraded  and/or  slowly  destabilizing  performance.  This  also  supports 
the  fact  that  the  LPV  controller  has  to  focus  the  majority  of  its  efforts  on  stability, 
as  was  found  earlier. 
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The  short  period  LTI  controller  exhibits  similar  behavior  when  simulated  with 
the  full  longitudinal  simulation  model;  however,  it  was  designed  using  the  short 
period  model.  Figure  5.25  shows  the  response  it  would  have  at  the  center  and  at 
corners  of  the  LPV  design  envelope  when  tested  with  its  design  plant.  As  expected,  at 
the  center  of  the  envelope  (its  design  point)  the  controller  has  excellent  performance; 
elsewhere,  it  appears  to  be  stable  over  a  large  portion  of  the  envelope,  though  at 
the  cost  of  some  tracking  error.  Nevertheless,  Figure  5.25  supports  the  fact  that, 
unlike  the  full  longitudinal  LPV  controller,  the  short  period  LPV  controller  does 
not  have  to  spend  as  much  effort  to  guarantee  robust  stability,  and  can  thus  devote 
more  effort  to  achieving  perfomance.  Note  also  that  several  LTI//i  controllers  could 
give  much  better  performance  when  simulated  with  the  short  period  plant;  however, 
this  controller  was  chosen  since  it  gave  better  performance  when  tested  with  the  full 
longitudinal  simulation  model. 

Finally,  this  short  period  LTI  controller  was  used  in  a  dynamic  simulation  to 
compare  its  performance  to  that  of  the  short  period  LPV  controller.  Its  response  to 
the  doublet  pitch-rate  command  is  shown  in  Figures  5.26  to  5.29,  while  that  of  the 
LPV  controller  is  given  in  Figures  5.11  to  5.14.  It  does  surprisingly  well,  although  it 
does  exhibit  some  diverging  behavior  after  extended  periods  of  time,  slightly  more 
overshoot,  and  tracking  is  not  quite  as  good.  While  the  LPV  controller  definitely 
performs  better,  this  does  illustrate  that  a  well-designed  LTI  controller  may  perform 
nearly  or  just  as  well  as  an  LPV  controller.  However,  since  it  does  directly  design  a 
varying  controller,  the  LPV  design  method  should,  in  general,  perform  better.  The 
extent  to  which  it  does,  in  all  likelihood,  depends  largely  on  the  extent  the  plant 
varies  and  the  rate  at  which  its  parameters  change.  For  linear  systems,  these  are 
aspects  which  only  this  or  similar  LPV  approaches  to  gain-scheduling  can  hope  to 
address. 
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Figure  5.24  Step  responses  using  full  longitudinal  LTI  controller 
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Figure  5.26  Dynamic  response  (using  LTI  controller)  to  a  positive  doublet  pitch- 
rate  command:  qc,  q 
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Figure  5.29  Dynamic  response  (using  LTI  controller)  to  a  positive  doublet  pitch- 
rate  command:  a,  6 
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VI.  Conclusions  and  Recommendations 


The  objectives  of  this  thesis,  set  forth  in  Chapter  I,  were  accomplished.  The 
LMI  formulation  of  a  new  gain-scheduling  approach  for  LPV  systems  developed  by 
Packard,  Apkarian,  and  Cabinet  was  presented,  along  with  the  supporting  theory 
and  an  extension  to  address  LPV  systems  with  modelled  uncertainty.  The  method 
was  tested  in  a  realistic  aircraft  control  problem  in  order  to  design  two  parameter- 
varying  pitch-rate  controllers,  one  for  the  short  period  aircraft  model  and  the  second 
for  the  full  longitudinal  model.  The  resulting  short  period  LPV  controller  gives  the 
desired  robust  performance  and  smoothly  schedules  itself  based  on  measured  values 
of  the  operating  parameters.  In  the  case  of  the  full  longitudinal  design,  the  resulting 
LPV  controller  is  robustly  stabilizing,  but  does  not  provide  acceptable  performance. 
This  emphasizes  the  conservatism  of  the  approach  which  greatly  increases  when  for¬ 
mulating  large  problems  with  many  parameter-dependent  elements,  such  as  the  full 
longitudinal  design  problem.  Nevertheless,  in  the  light  of  the  simulation  results,  the 
short  period  LPV  controller  appears  to  be  sufficient  to  provide  the  desired  responses 
in  the  design  envelope,  although  full  non-linear  aircraft  models  should  be  tested  to 
confirm  this. 

6. 1  Conclusions 

Based  on  the  presented  theory  and  the  subsequent  implementations,  some  con¬ 
clusions  can  readily  be  drawn. 

First  and  foremost  is  the  fact  that  this  new  method  provides  a  relatively 
straightforward  and  numerically  tractable  way  to  design  a  practical  gain-scheduled 
controller.  As  opposed  to  more  classical  methods,  no  point  designs  are  required  and 
no  schedule  must  be  determined;  instead,  the  control  law  is  directly  designed  by 
solving  a  convex  optimization  problem.  In  addition,  unlike  the  other  similar  ap- 
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proaches  mentioned  in  Chapter  I  where  an  infinite  number  of  LMI’s  must  be  solved 
[BP94],[WYPB94],  this  method  only  requires  solving  4  LMPs. 

However,  as  has  been  explained,  depending  on  the  size  of  the  design  problem, 
even  solving  these  four  LMPs  can  prove  to  be  a  difficult  task,  since  the  LMI  solvers 
are  extremely  computationally  demanding.  The  size  of  the  problem  can  be  reduced 
by  limiting  the  order  of  the  parameters  on  which  the  elements  of  the  state-space 
matricies  depend,  but  this  comes  at  the  expense  of  additional  model  uncertainties. 

In  addition,  using  this  method,  it  is  unlikely  that  a  single  LPV  controller  can 
be  found  to  gain-schedule  an  entire  aircraft  operating  envelope.  This  is  undoubtedly 
due  to  the  excessive  conservatism  in  the  approach  which  stems  from  several  sources. 
While  it  simplifies  the  formulation  and  solvability  of  the  problem,  the  Small  Gain 
Theorem  forces  the  controller  to  handle  complex,  usually  non-existing  values  of  the 
parameters  in  addition  to  the  real  ones.  Also,  the  problem  formulation  allows  for 
infinitely  fast  parameter  variation  rates.  Finally,  if  ^-synthesis  is  used,  the  fact 
that  only  complex  full  blocks  can  be  used  to  represent  the  parameter  blocks  further 
increases  the  conservatism  of  the  resulting  controller  designs.  The  result  of  this 
compounded  conservatism  is  undoubtedly  degraded  robust  performance,  and  results 
m  guaranteed  robust  performance  in  a  smaller  envelope  than  is  likely  achievable  in 
practice.  For  this  reason,  testing  is  recommended  beyond  the  guaranteed  envelope 
to  further  stretch  the  envelope  over  which  the  LPV  controller  is  effective. 

The  last  paragraph  notwithstanding,  the  D-K-D  procedure  described  in  the 
thesis  is  strongly  recommended  for  all  gain-scheduling  design  problems,  even  when 
additional  uncertainty  is  not  modelled  in  the  design,  because  the  scales  resulting  from 
the  gain-scheduling  part  of  the  problem  only  optimize  the  design  on  the  input /output 
channels  corresponding  to  the  varying  parameter  block,  whereas  the  scales  derived 
in  //-synthesis  are  meant  to  optimize  the  design  on  all  the  input/output  channels. 
In  any  case,  the  controller  for  the  scaled  H^o  problem  (the  original  gain-scheduling 
problem)  is  still  available  from  the  first  D-K-D  iteration. 
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It  was  also  demonstrated,  in  the  case  of  the  short  period  design,  that  an  LPV 
controller  could  guarantee  robust  stability  and  nominal  performance  over  a  much 
wider  envelope  than  any  optimal  LTI  controller.  This  was  especially  significant  for 
a  more  complex  plant  like  that  of  the  full  longitudinal  design,  for  which  the  LTI 
norm  for  a  robustly  stabilizing  LPV  envelope  reached  a  whopping  157.41.  Having 
said  that,  the  tradeoff  between  performance  and  stability  is  still  present  and,  as 
the  full  longitudinal  design  demonstrated,  may  compete  to  the  point  that  both  re¬ 
quirements  cannot  be  met  simultaneously.  This  then  requires  either  reducing  the 
design  envelope,  and  possibly  simplifying  the  problem,  or  changing  the  performance 
requirements. 

Specifically,  if  further  full  longitudinal  designs  are  considered,  it  may  be  neces¬ 
sary  to  include  a  more  complex,  possibly  time- varying,  ideal  reference  model;  alter¬ 
natively,  the  model-matching  design  approach  could  be  replaced  by  a  design  model 
in  which  performance  is  maximized  by  weighting  the  sensitivity.  This  last  approach 
would  avoid  the  need  for  a  model  and  simply  provide  the  best  performance  avail¬ 
able.  Another  possibility  would  be  to  minimize  the  number  of  parameter-dependent 
stability  derivatives  by  using  a  constant  value  for  those  that  have  the  least  effect  on 
the  overall  aircraft  model.  Whatever  the  design  methodology  adopted,  the  designer 
should  first  develop  the  best  LTI  design  possible.  At  best,  this  may  eliminate  the 
need  to  design  an  LPV  controller  if  performance  is  acceptable.  At  worst,  it  will 
guide  the  proper  choices  of  weights  required  to  optimize  the  problem.  In  addition, 
the  norm  of  the  LTI  nominal  performance  problem  should  be  tested  at  the  center  of 
the  robustly  stabilizable  design  envelope,  as  an  indicator  of  whether  or  not  robust 
performance  will  likely  be  achievable. 

Finally,  it  has  been  shown  that,  for  at  least  some  control  problems,  LTI  con¬ 
trollers  may  perform  as  well  or  nearly  as  well  as  LPV  controllers  designed  by  this 
method.  This  undoubtedly  depends  upon  the  complexity  of  the  design  problem  and 
the  variability  of  the  plant  and  its  operating  conditions.  Since  the  LPV  controllers 
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adapt  to  the  changing  operating  conditions,  their  performance  should  in  general  sur¬ 
pass  that  of  equivalent  LTI  controllers  when  the  plant  varies  significantly  throughout 
its  operating  envelope. 

6.8  Recommendations 

While  the  results  were  encouraging,  some  of  the  shortcomings  of  this  gain¬ 
scheduling  method  highlight  the  need  for  further  research  in  some  areas. 

Above  all,  there  is  the  need  to  reduce  the  amount  of  conservatism  inherent 
in  the  present  method.  This  can  be  done  by  reducing  the  effect  of  or,  preferably, 
entirely  eliminating  one  or  more  of  the  three  sources  of  conservatism  mentioned 
in  the  last  section.  Becker’s  method  [BP94],  which  considers  only  real  values  of 
the  time- varying  parameters,  should  be  investigated  and  compared  to  the  results 
obtained  using  the  current  method.  Since  it  involves  discretizing  an  infinite  number 
of  constraints,  a  small  design  problem  such  as  the  short  period  model  is  strongly 
recommended  since  the  LMI  solver  is  so  computationally  demanding.  The  method 
of  Wu  [WYPB94],  which  would  further  eliminate  infinitely  fast  parameter  variations, 
is  only  recommended  if  a  more  numerically  tractable  variation  of  the  method  can  be 
presented;  as  is,  even  a  small  problem  like  the  short  period  model  would  be  extremely 
demanding.  However,  a  possible  alternative  that  should  be  investigated  would  be  to 
include  parameter  rates  as  part  of  the  time- varying  parameter  block  by  incorporating 
them  in  the  curve-fits  for  the  elements  of  the  state-space  matricies.  This  would 
prevent  infinitely  fast  rates  of  change,  although  it  would  still  allow  complex  values 
for  them.  Generally,  though,  it  should  help  decrease  some  of  the  conservatism. 

The  advent  of  a  practical  method  for  /^-synthesis  of  real  scalar  uncertainties 
would  eliminate  the  third  source  of  conservatism.  In  fact,  by  combining  the  pa¬ 
rameter  block  for  the  LPV  controller  and  that  of  the  LPV  plant  and  by  grouping 
like  parameters  to  form  real  repeated  scalar  blocks,  such  real  /x-synthesis  techniques 
should  entirely  eliminate  the  need  for  the  LMI  Li-scales,  since  these  would  effectively 
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already  be  produced  by  the  /i-synthesis  iteration.  Certainly,  when  real  //-synthesis 
becomes  practical,  gain-scheduling  should  thus  be  revisited. 

An  observation  can  also  be  made  that,  while  the  more  classical  methods  must 
determine  control  laws  for  the  controller  or  its  coefficients  (depending  on  the  method), 
the  LPV  method  discussed  in  this  thesis  requires  determining  a  curve-fit  for  many 
elements  of  the  state-space  matricies.  However,  the  LPV  method  is  able  to  guarantee 
robust  performance  over  its  entire  design  envelope,  whereas  the  others  cannot.  This 
also  emphasizes  the  importance  of  maintaining  relatively  good  accuracy  in  mod¬ 
elling  the  LPV  plant;  otherwise,  robust  performance  will  be  meaningless.  This  is 
an  issue  that  warrants  further  investigation;  specifically,  to  examine  the  relationship 
between  curve-fitting  the  elements  of  the  LPV  state-space  matricies,  the  resulting 
plant  inaccuracies,  and  the  effect  on  actual  controller  robust  performance. 

Finally,  further  controller  designs  should  also  be  attempted  in  discrete-time 
and  account  for  time  delays  and  errors  in  the  parameter  measurements  which  would 
occur  in  an  actual  aircraft  implementation.  The  discrete-time  problem  formulation  is 
already  provided  in  [AG95].  Similarly,  incorporating  the  advantages  of  other  norms 
such  as  H2  and  /i,  to  reduce  the  effects  of  noise  and  add  limits  on  control  deflection 
or  error  magnitude  should  also  be  considered.  An  LMI  approach  for  a  general  mixed 
H2IH00  problem  has  already  been  presented  in  [CW94];  similarly,  a  mixed  H2IH00 
approach  has  been  applied  to  LPV  systems  in  [Sch95].  These  methods  could  be 
further  extended  to  the  gain-scheduled  LPV  Hoo  control  method  discussed  in  this 
thesis. 
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Appendix  A.  F-18  Design  Flight  Conditions 

The  following  table  lists  the  trimmed  flight  conditions  that  were  used  in  the 
design  and  simulation  of  the  longitudinal  gain-scheduled  controller.  The  data  point 
numbers  correspond  to  those  given  in  all  the  figures  of  Appendix  B,  in  order  to  easily 
cross-reference  them  with  the  stability  derivative  curve-fits. 


Table  A.l  Longitudinal  Design  Flight  Conditions 


Data 

Point 

# 

Mach 

Number 

Altitude 

(ft) 

Dynamic 

Pressure 

(psf) 

a 

(deg) 

1 

0.3 

26000 

47.4 

25.2 

2 

0.5 

40000 

68.5 

16.8 

3 

0.6 

30000 

158.4 

5.2 

4 

0.4 

6000 

189.9 

6.0 

5 

BESsIsH 

2.6 

6 

1.9 

7 

1.6 

8 

0.8 

10000 

1.7 

9 

0.8 

5000 

1.5 

■n 

0.9 

10000 

825.2 

1.4 

■B 

0.85 

5000 

890.8 

1.4 

mm 

0.9 

5000 

998.7 

1.3 

■B 

15000 

75.0 

15.0 

14 

12000 

236.0 

4.8 

mm 

0.5 

20000 

6.5 

■B 

0.6 

5000 

2.6 

17 

0.8 

25000 

3.0 

18 

0.8 

35000 

125.0 

4.5 
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Appendix  B.  F-18  Stability  Derivative  Curve  Fits 

This  appendix  illustrates  the  curve-fits  used  to  establish  relationships  between 
the  stability  derivatives  and  the  varying  parameters  M  and  h  for  the  simulation 
model,  the  short  period  design  plant,  and  the  full  longitudinal  design  plant.  For 
all  the  figures,  the  data  point  numbers  correspond  to  those  given  in  Table  A.l  of 
Appendix  A,  in  order  to  easily  cross-reference  the  actual  trimmed  flight  conditions. 

B.l  Simulation  Model  Curve  Fits 

The  following  figures  compare  the  trimmed  data  points  to  the  polynomial  curve 
fits  used  for  the  simulation  model  stability  derivatives.  The  polynomials  are  as  high 
as  order  and  cover  the  entire  available  flight  envelope  such  that 

M  6  [0.3, 0.95]  and  h  E  [5000, 40000]  ft.  (B.l) 


Figure  B.l  Simulation  model  curve  fit  of  Xu 
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Figure  B.7  Simulation  model  curve  fit  of  Z, 


Zdelta 


2 


4 


6 


B.2  Short  Period  Design  Model  Curve  Fits 

The  following  figures  compare  the  trimmed  data  points  to  the  polynomial  curve 
fits  used  for  the  short  period  design  model  stability  derivatives.  The  polynomials  are 
no  higher  than  order  M'^h  and  the  design  envelope  is 

M  6  [0.5, 0.95]  and  h  6  [5000, 30000]  ft.  (B.2) 


Figure  B.16  Short  period  design  model  curve  fit  of  Za 
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Figure  B.19  Short  period  design  model  curve  fit  of 
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Figure  B.20  Short  period  design  model  curve  fit  of  Mq 


Mdelta 


B,3  Full  Longitudinal  Design  Model  Curve  Fits 


The  following  figures  compare  the  trimmed  data  points  to  the  polynomial  curve 
fits  used  for  the  full  longitudinal  model  stability  derivatives.  The  polynomials  are 
no  higher  than  order  M'^h  and  the  design  envelope  is 

M  6  [0.5, 0.95]  and  h  6  [5000, 30000]  ft.  (B.3) 


Figure  B.22  Full  longitudinal  design  model  curve  fit  of  Xu 


Appendix  C.  Implementation  Programs 


This  appendix  lists  the  main  implementation  programs,  along  with  the  main 
subprograms.  They  are  included  here  for  informational  purposes  only;  the  author 
should  be  contacted  for  access  to  the  original  set  of  programs  (which  are  always 
evolving!). 

C.l  fl81ongdat.ni;  F-18  Trimmed  Flight  Condition  Data 


•/•/•/•/ 0/ i/Vl/ •/•/•/•/O/VVO/ •/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/«/ VVVVV'i 


'/,  flSlongdat.m 

y,  F18  data  for  linear  longitudinal  models 


y,  m3h26 

al=[-8.6469e-4  3.4079e-2  -2.242e0  -5.1786e-l; 
-3.6574e-2  -2.2959e-l  9.9312e-l  3.9494e-2; 
-1.6123e-2  2.4361e-2  -2.0464e-l  2.0935e-3; 

0  0  9.998e-l  0] ; 

bl= [-5 . 2 165e-2 ; -4 . 0340e-2 ; - 1 . 7302 ; 0] ; 
y,  m5h40 

a2=[-2.4351e-3  6.204e-3  -2.4285e0  -5.4746e-l; 
-1.5514e-2  -2.4232e-l  9.9641e-l  1.9699e-2; 
-1.9554e-2  -2.3421e0  -1.7374e-l  -6.605e-4; 

0  0  9.9998e-l  0] ; 

b2= [9 . 3859e-3 ; -4 . 1597e-2 ; -2 . 5953 ; 0] ; 
y.  m6h30 

b3= [-9 . 1427e-3 ; -9 . 2769e-2 ; -6 . 5735 ; 0] ; 
a3=[-4.0006e-3  -6.6825e-2  -1.2218e0  -5.5765e-l; 
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-1.1737e-2  -5.0877e-l  9.9395e-l  2.1962e-3; 
-1.7155e-2  -1.1314eO  -2.8045e-l  -1.0661e-2; 

0  0  10]; 

•/,  m4h6 

b4= [2 . 7435e-2 ; -1 . 5077e- 1 ; -7 . 9255 ; 0] ; 

a4=[-3.4565e-3  -2.3318e-2  -7.8631e-l  -3.7581e-l 
-2.0661e-2  -8.0179e-l  9.8467e-l  1.6962e-2; 
-7.9341e-3  -1.5211e0  -5.944e-l  1.9196e-l; 

0  0  10]; 

•/.  m7hl4 

b5= [-1 .4391e-l ; -1 . 9397e-l ; -1 . 9292el ; 0] ; 

a5=[- 1.29476-2  5.6579e-l  -5.8325e-l  -5.2547e-l; 
-7.9393e-3  -1.1751e0  9.8710e-l  8.2546e-3; 
5.9347e-3  -8.4584e0  -8.7762e-l  3.4102e-3; 

0  0  10]; 

•/.  m8hl2 

b6= [-1.72016-1 ; -2.31576-1 ; -2 . 647961 ; 0] ; 

a6=[-l. 75116-2  1.558660  -4.74176-1  -5.65866-1; 
-6.02656-3  -1.562460  9.8619e-l  8.79446-3; 
-1.56836-2  -1.493961  -1.132160  8.26646-3; 

0  0  10]; 

•/.  m95h20 

b7= [1 .4969e-l ; -1 . 86726-1 ; -2 . 721661 ; 0] ; 

a7= [-7. 55066-2  6.48046-1  -4.6379e-l  -5.6189e-l; 
-7.82766-3  -1.905460  9.8956-1  1.50976-2; 
-3.19256-1  -3.388561  -9.87166-1  1.1269e-l; 

0  0  10]; 

•/.  m8hl0 
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b8= [-2 . 9175e-l ; -2 .4491e-l ; -2 . 8335el ; 0] ; 

a8=[-2.1964e-2  8.823e-l  -4.5306e-l  -5.7631e-l; 
-5.7642e-3  -1.6750e0  9.9531e-l  1.1595e-2; 

-9.8837e-3  -1.6158el  -1.2120e0  -1.8555e-2; 

0  0  1  0]  ; 

•/.  m8h5 

b9=  [-1 . 9651e-l ;-2 .8517e-l ;-3 .3445el ; 0] ; 

a9=[-2.4285e-2  1.8452e0  -3.8628e-l  -5.6327e-l; 
-5.1311e-3  -1.9937e0  9.8277e-l  1.2773e-2; 
5.6171e-3  -1.9439el  -1.4272e0  6.5051e-3; 

0  0  10]; 

•/,  m9hl0 

bl0= [-4 . 4621e-l ; -2 . 7569e-l ; -3 . 7363el ; 0] ; 

al0=[-4.4422e-2  8.6444e-l  -4.2435e-l  -4.0672e-l 
-9.5576e-3  -2.4524e0  9.8565e-l  1.7780e-2; 
-2.4665e-l  -3.8613el  -1.3401e0  1.9235e-l; 

0  0  10]; 

•/.  ni85h5 

bl 1= [-3 . 9377e- 1 ; -3 . 0 120e-l ; -3 . 8430el ; 0] ; 

all=[-3.5431e-2  1.1270e0  -3.8403e-l  -5.7154e-l; 
-3.6673e-3  -2.3281e0  9.8312e-l  1.1098e-2; 
1.8523e-2  -3.0436el  -1.4930e0  3.5754e-2; 

0  0  10]; 

'/,  m9h5 

al2=[-.0516  1.3906  -0.3817  -0.5485; 

-0.0094  -2.9105  0.9835  0.0141; 

-0.2600  -46.4719  -1.5527  -0.0948; 

0  0  10]; 

bl2= [-0.4662; -0.3161; -43. 6515 ;0] ; 
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•/.m3hl5 


al3  =[-9.1762e-04  -4.2476e-02  -1.4162e+00  -5.4260e-01 

-3.5627e-02  -3.3711e-01  9.8974e-01  1.5817e-05 

2.2474e-03  -8.2586e-01  -2.9308e-01  -2.9711e-05 

0  0  9.9998e-01  0]; 


bl3  =[-3.9377e-02 
-6.8546e-02 
-2.6358e+00 
0]; 

•/,m5hl2 

al4  =[-1.6397e-03  2.7965e-01  -7.7173e-01  -5.5977e-01 

-1.3957e-02  -8.2465e-01  9.8680e-01  8.2671e-05 

-2.0111e-02  -1.8828e+00  -6.0418e-01  1.5988e-04 

0  0  l.OOOOe+00  0]; 


bl4  =[3.5500e-02 
-1.5267e-01 
-9.7970e+00 
0]; 


•/,m5h20 

al5  =[-3.7742e-06  1.3294e-02  -l.OlOle+00  -5.5813e-01 

-1.4848e-02  -6.1826e-01  9.9076e-01  4.7958e-06 

-1.9803e-02  -1.3457e+00  -4.0138e-01  -6.1538e-07 

0  0  l.OOOOe+00  0]; 


bl5  =[1.5361e-02 
-1.1430e-01 
-7.0813e+00 
0]; 
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'/,m6h5 


al6  =[-7.3877e-03  5.6551e-01  -5.2056e-01  -5.6125e-01 

-9.6253e-03  -1.2925e+00  9.8260e-01  8.9218e-05 

4.8902e-03  -8.5064e+00  -9.7267e-01  1.1681e-05 

0  0  l.OOOOe+00  0]; 


bl6  =[-1.4960e-01 
-2.2700e-01 
-1.9781e+01 
0]; 


•/,m8h25 

ai7  =[-4.1811g-03  7.9322e-01  -7.2318e-01  -5.6097e-01 

-7.0760e-03  -9.3598e-01  9.9135e-01  6.3784e-06 

-3.2723e-02  -7.2386e+00  -7.0232e-01  2.1230e-04 

0  0  l.OOOOe+00  0]; 


bl7  =[5.6906e-02 
-1.4529e-01 
-1.5594e+01 
0]; 


•/.m8h35 

al8  =[-1.3032e-03  2.8386e-01  -1.0511e+00  -5.5998e-01 

-8.0363e-03  -6.2303e-01  9.9446e-01  1.2041e-05 

-4.4622e-02  -4.1277e+00  -4.5340e-01  1.6097e-04 

0  0  l.OOOOe+00  0]; 


bl8  =[l.4165e-02 
-9.7770e-02 
-9.9299e+00 
0]; 
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'/,  M,h,qbar(dyn  press), aoa  for  all  12  points 


Mhqa=[.3  26  47.4  25.2  ;  .5  40  68.5  16.8 
.6  30  158.4  5.2  ;  .4  6  189.9  6 
.7  14  426.4  2.6  ;  .8  12  603  1.9 
.95  20  614.4  1.6  ;  .8  10  652  1.7 
.8  5  789.1  1.5  ;  .9  10  825.2  1.4 
.85  5  890.8  1.4  ;  .9  5  998.7  1.3]; 


Mhqa2=[.3  15  75  15.0 


5 

5  20 

12 

170 

236 

6.5 

4.8 

6 

5 

444 

2.6 

8 

25 

198 

3.0 

8 

35 

125 

4.5]  ; 

Mhqa=[Mhqa;Mhqa2] ; 
data=size(Mhqa, 1) ; 


'/,plot(Mhqa(:  ,1)  ,Mhqa(:  ,2)  , 'x')  ‘/.plot  M  vs  h 
'/,axis([.2  1  0  50]); 

'/.xlabel  ( '  Mach  number  ’ )  ; 

'/.ylabelC 'altitude  (1000  ft)'); 


y,  end  of  fl81ongdat.m 


e/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ Q/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ S/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ P/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ 0/ a/ 0/ Oi 


C.2  fl8poly.m;  Curve-fit  Trimmed  Data  to  Polynomials 


•/•/•/•/•/•/•/•/•/ 0/0/ O/O/I/I/O/ •/•/I/O/ •/•/•/o/o/vi/o/ 0/0/ •/0/0/vwo/o/vVO/VVVO/o/VVO/VV‘J 


y,  fl8poly.m 

y,  converts  full  longitudinal  data  to  polynomial  curve  fit 
f 181ongdat ; 
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'/,  choose  design  envelope 
Mmin= . 5 ; 

Mmax=.95 ; 
hmin=5 ; 
hmax=30 ; 

*/,  build  Adat  for  all  pts  (combines  a,b) 

Adat=[]  ; 
pts=0; 

M=[]; 
h=[]; 
dindx=[]  ; 

for  k=l:data, 

Mdat=Mhqa(k,l) ; 
hdat=Mhqa(k,2) ; 

if  (Mdat>=Mmin  &  Mdat<=Mmax  &  hdat<=hmax:  &  hdat>=hmin) 
pts=pts+l ; 

M=[M;Mdat] ; 
h=[h;hdat] ; 
dindx=[dindx;k] ; 
eval(C'a=a'  int2str(k) 
eval(['b=b'  int2str(k) 

Adat=[Adat;  a(l,:)  b(l)  a(2,:)  b(2)  a(3,:)  b(3)] ; 
end 
end 

clear  a  b 

qbar=Mhqa(: ,3) ; 
qi=diag(minv(diag(qbar))) ; 
qi2=qi.*qi; 
q2=qbar . *qbar ; 

•/.h=1000*Mhqa(:  ,2)  ; 

h=1000*h; 

h2=h.*h; 

•/.M=Mhqa(:,l); 

Mi=diag(minv(diag(M))) ; 

Mi2=Mi.*Mi; 

Mh=M.*h; 

M2=M.*M; 
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M2h=M2.*h; 

Mh2=M.*h2; 

M2h2=M2.*h2; 

M3=M2.*M; 

M3h=diag(M3*h'); 

M3h2=diagCM3*h2') ; 


'/,  determine  relationships  between  Adat  and  M,h 

*!,  names  of  stab,  derivatives  that  depend  on  M,h 

names* ' Xu  Xa  Xq  Xt  Xde  Zu  Za  Zq  Zt  Zde  Mu  Ma  Mq  Mt  Mde ' ; 

[num , names , namelen ,maxlen] =namstuf f (names) ; 

y,  most  accurate  for  generating  only  10  terms 

order* [2  1;2  1;2  1;2  1;2  1;1  1;1  1;1  1;1  1;1  1;2  1;2  1;2  1;2  1;1  1]; 
y  Dr  R. ' s  eyeball 

y.order*[2  1;3  2;1  1;3  2;3  2;2  1;1  1;0  0;2  1;1  1;3  2;2  1;1  1;3  2;1  1]  ; 
y,  <*  25'/.  error 

'/.order=[3  2;3  2;2  1;2  2;3  2;3  2;1  1;0  0;2  2;2  1;3  2;3  2;1  2;3  2;1  1]  ; 
'/,  <*  10'/,  error 

•/,order=[3  2;3  2;2  2;2  2;3  2;3  2;2  1;0  0;3  2;2  1;3  2;3  2;1  2;3  2;2  1]  ; 
'/.  <*  1'/.  error 

'/.order=[3  2;3  2;2  2;3  2;3  2;3  2;3  1;0  0;3  2;3  1;3  2;3  2;2  2;3  2;3  1]  ; 
y,  <=  ly,  error  and  delta  s.v.  <  -20dB 

'/.order=[3  2;3  2;2  2;3  2;3  2;3  2;3  2;1  1;3  2;3  1;3  2;3  2;2  2, -3  2;3  2]; 
y,  MS  order 

border* [2*ones (15,1)  ones(15,l)] ; 

y,  superduper  —  used  for  simulation  model 
sim_order*[3*ones(15,l)  2*ones(15,l)] ; 
yorder=[3*ones(15,l)  2*ones(l5,l)] ; 

•/.  q 

border* Cones (15,1)  zeros(l5,l)] ; 
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'/,  constaint  nominal  plant 
'/,order=zeros(15,2) ; 

maixfit=max( [order ( :  ,1)  ;order(:  ,2)])  ; 

aaOO=[ones(pts,l)] ; 
aall=[aa00  h  M  Mh] ; 
aa21=[aall  M2  M2h] ; 
aa22=[aa00  h  h2  M  Mh  Mh2] ; 
aal2=[aa00  h  h2  M  Mh  Mh2] ; 
aa22=[aal2  M2  M2h  M2h2] ; 
aa31=[aa21  M3  M3h] ; 
aa32=[aa22  M3  M3h  M3h2] ; 

'/,aal0=[aa00  qbar]  ; 

for  k=l:num; 

□M=order(k, 1) ; 

Dh=order(k,2) ; 

eval(['aa=aa'  int2str(0M)  int2str(0h) 

name=names (k , 1 : namelen (k) ) ; 

bb=Adat ( : , k) ; 

cc=aa' *aa; 

dd=aa' *bb; 

xx=cc\dd; 

eval ( [ ’ P0LY_ ‘  name  ’  =xx ; ' ] ) ; 

err=[3  ; 

for  i=l:pts, 

er=(bb(i)-aa(i, : )*xx)/bb(i) ; 
err= [err ; er] ; 
end 

err=max(abs(err)) ; 
eval(['err_’  name  '=100*err'] ) ; 
end 

'/,  plots 
plt=0; 

if  plt==l 

for  k=l:num, 
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OM=order(k,l) ; 

0h=order(k,2) ; 

aval ( [ ' aa=aa '  int2str(DM)  int2str(0h) 
name=neLmes(k,l:namelen(k)) ; 
t=l:pts; 
figure (k) 

plot (dindx(t) ,Adat(: ,k)  ,  'x') 
legend ('data  point'); 
hold; 

aval ( ['plot (dindx(t) ,aa*P0LY_'  name  ')']); 
y,  legend  ('curve  fit'); 
hold; 

y,  title  (name); 

xlabel('data  point  #'); 
ylabel(name)  ; 
end 


end 

y,  end  of  flSpoly.m 


C.3  flSsppoly.m;  Curve-fit  Trimmed  Data  to  Polynomials 


VVVVVVVVVVVVVVVVVVVVVVVV*/VVVVVVO/o/V*/*/*/o/c/o/i/»/«/o, 


y,  flSsppoly.m 

%  converts  longitudinal  data  to  short  period  polynomial  curve  fit 
f ISlongdat; 

7,  choose  design  envelope 
Mmin=.4; 

Mmax=.95; 
hmin=5 ; 
hmax=30 ; 

y  build  Adat  for  all  pts  (combines  a,b) 
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Adat=[]  ; 
pts=0; 

M=[]; 
h=[]; 
dindx=  [] ; 

for  k=l :data, 

Mdat=Mhqa(k, 1) ; 
hdat=Mhqa(k,2) ; 

if  (Mdat>=Mmin  &  Mdat<=Mmax  &  hdat<=hm2LX  &  hdat>=hmin) 
pts=pts+l ; 

M=[M;Mdat] ; 
h= [h;hdat] ; 
dindx= [dindx ;k] ; 
eval(C'a=a'  int2str(k) 
eval(['b=b'  int2str(k) 

Adat=[Adat;  a(2,2:3)  b(2)  a(3,2:3)  b(3)] ; 
end 
end 

clear  a  b 

y,  determine  relationships  between  Adat  and  M,h 

*/,  names  of  stab,  derivatives  that  depend  on  M,h 
names='Xu  Xa  Xq  Xt  Xde  Zu  Za  Zq  Zt  Zde  Mu  Ma  Mq  Mt  Mde^ ; 

[num ,  names ,  namelen , maxilen]  =namstuf f  (names )  ; 

qbar=Mhqa( : , 3)  ; 
qi=diag(minv(diag(qbar))) ; 
qi2=diag(qi*qi') ; 
q2=diag(qbar*qbar' ) ; 
y,h=1000*Mhqa(:  ,2) ; 
h=1000*h; 

h=(h-22500) ./35000; 

h2=diag(h*h'); 

y.M=Mhqa(:.l); 

M=(M-.625) ./.65; 

Mi=diag(minv(diag(M) ))  ; 

Mi2=diag(Mi*Mi') ; 
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Mh=diag(M*h' ) ; 

M2=diag(M*M0; 

M2h=diag(M2*hO; 

Mh2=diag(M*h20; 

M2h2=diag(M2*h2'); 

M3=diag(M2*MO; 

M3h=diag(M3*hO; 

M3h2=diag(M3*h2') ; 


'/,  determine  relationships  between  Adat  and  M,h 

'/,  names  of  stab,  derivatives  that  depend  on  M,h 
names='Zalpha  Zq  Zdelta  Malpha  Mq  Mdelta' ; 

[num , neimes ,namelen,maxlen]  =namstuff  (neimes)  ; 

'/,  most  accurate  for  generating  only  6  terms 
order=[l  1;1  1;1  1;2  1;2  1;1  1]; 

y,  <=  25y,  error  —  10  term  theta  block 
y.order=[l  1;0  0;2  1;3  2;1  2;1  1]; 

y,  <=  lOy,  error  —  10  terms 
y.order=  [2  1;0  0;2  1;3  2;1  2;2  1]; 

'/i  <=  ly  error  —  11  terms 
y,order=[3  1;0  0;3  1;3  2;2  2;3  1]; 

y,  <=  ly  error  and  delta  s.v.  <  -20dB  —  12  terms 

y.order=[3  2;1  1;3  2;3  2;3  2;3  2]; 

y,  7  terms 

y,order=[2*ones(6,l)  ones (6,1)]  ; 
y,  8  terms 

y.order=[l  1;0  0;1  1;3  2;1  1;2  1]; 

y,  modified  MS  —  6  term  theta 
y,order=[l  1;0  0;1  1;2  1;1  1;1  1]; 

y,  superduper  —  12  terms 
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'/,or<ier=[3*ones(6,l)  2*ones(6,l)]  ; 
y,order=  [3  2;1  1;3  2;3  2;3  2;3  2]; 

•/.  q 

'/,ordGr=  [ones (6,1)  ZGros(6,l)]; 

'/,  constant  nominal  plEint 
y,order=zeros(6,2) ; 

maxfit=max([order(: ,1) ;order(: ,2)]) ; 

aaOO=[ones(pts,l)] ; 
aall=[aa00  h  M  Mh] ; 
aa21=[aall  M2  M2h] ; 
aa22=[aa00  h  h2  M  Mh  Mh2] ; 
aal2=[aa00  h  h2  M  Mh  Mh2] ; 
aa22=[aal2  M2  M2h  M2h2] ; 
aa31=[aa21  M3  M3h] ; 
aa32=[aa22  M3  M3h  M3h2] ; 

y,aal0=[aa00  qbar]  ; 

for  k=l:num; 

QM=order(k,l) ; 

0h=order(k,2) ; 

eval( ['aa=aa'  int2str(0M)  int2str(0h) 

naLme=names  (k ,  1 :  namelen  (k)  )  ; 

bb=Adat ( : , k) ; 

cc=aa' *aa; 

dd=aa'*bb; 

xx=cc\dd ; 

eval(['P0LY_'  name  '=xx;']); 

err=[]  ; 

for  i=l:pts, 

er=(bb(i)-aa(i, : )*xx)/bb(i) ; 
err= [err ; er] ; 
end 

err=max(abs(err)) ; 
eval(['err_'  name  '=100*err']) ; 
end 
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plots 

plt=l; 

if  plt==l 

for  k=l :num, 

0M=order(k,l) ; 

0h=order(k,2) ; 

eval ( [ ’ aa=aa '  int2str(0M)  int2str(0h) 
naiiie=najaes(k,l:n2anelen(k))  ; 
t=l :pts; 
figure (k) 

plot  (dindx(t)  ,Adat  ( :  ,k)  , 'xO 
legend ('data  point'); 
hold; 

eval( ['plot(dindx(t) ,aa*P0LY_'  name  ')']); 
hold; 

*/,  title  (name); 

xlabel(' data  point  #'); 
ylabel (name) 
end 


end 

*/,  end  flSsppoly.m 


vir‘ 


C.4  f  18poly21ft  .m;  Convert  Polynomials  to  LFT 


*/,  f  18poly21ft .m 

•/. 

*/,  This  script  file  defines  the  parameter  normalization  factors, 
calculates  polynomial  coefficients  in  the  normalized  deltas, 
y,  and  constructs  LFTs  (in  deltas)  for  Za,  Zde,  Ma,  Mq,  Mde; 
y,  these  LFTs  are  combined  to  form  lftab=[Za  1  Zde;Ma  Mq  Mde] 
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Inputs:  names  =  a  string  containing  the  names  of  the  pletnt  parameters 
(Za,  Zq,  etc) 

POLY.xx  =  poly,  coef.s  for  the  parameters  given  in  names,  i.e. 

P0LY_Za,  POLY.Zq,  etc. 

maxfit  =  max.  order  of  the  poly's  defined  in  P0LY_xx 
order  =  matrix  defining  the  order  (in  M  and  h)  of  each  poly,  for 
each  of  the  derivatives  (i.e.  Za,  Zde,  ...) 

Each  row  of  order  coeresponds  to  one  of  the  derivatives, 
cind  the  are  listed  in  the  order  that  they  are  given  in 
names 

Outputs:  LFTs  for  each  of  the  derivatives  and  the  A  and  B  matrices. 

Written  by:  Capt  Martin  Breton,  AFIT  (for  full  longitudinal  +  schur  redux) 


y,  Based  on  a  short  period  program  written  by:  Capt  Mark  Spillman  and  Lawton  Lee 
y,  WL/FIGC-3  Bldg  146,  WPAFB 

y. 


y,  DEFINE  inputs 

n2unes='Xu  Xa  Xq  Xt  Xde  Zu  Za  Zq  Zt  Zde  Mu  Ma  Mq  Mt  Mde’ ; 

flag=order(8,l) ; 
if  flag  ==  0, 

order= [order (1 :7, :) ;order(9:15, :)] ; 

names= ‘ Xu  Xa  Xq  Xt  Xde  Zu  Za  Zt  Zde  Mu  Ma  Mq  Mt  Mde ’ ; 

LFT_Zq=P0LY_Zq; 

end 

tol=[]  ; 

Parameter  normalization  factors  and  shortcut  matrices  to  convert 
polynomial  coefficients  in  parm  to  coefficients  in  delta_parm 

y,  Design  and  spread  Machs 
Mbar= (Mmax+Mmin) / 2 ; 

Mtil=Mmax-Mbar ; 

y.  Design  cind  spread  altitudes 
hbar  =  500*(hmax+hmin) ; 
htil  =  1000*hmax-hbar; 

y,  Create  the  shortcut  matrices 
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Mmat=zeros(maxfit+l) ; 
hmat=zeros(niaLxfit+l)  ; 

K=zeros(mcLxf  it+1)  ;kvec=l ;K(1 , 1)=1 ; 
for  ii=l  :meLxfit+l; 

Mmat=Mmat  +  diag((Mtil.~(0:maxfit-ii+l))*(Mbar“(ii-l)) ; 
hmat=hmat  +  diag((htil.“(0:maxfit-ii+l))*(hbar"(ii-l)) ; 
if  ii<=meixfit 
kvec=conv(kvec, [1  -1]); 

K(l:ii+l,ii+l)=abs(kvec') ; 
end 
end 

Mmat=Mmat . *K ; 
hmat=hmat . *K ; 

y.mmy.y.y.  lft  contmction 

[niim,ncimes,namelen,maxlen]  =  neimstuff  (names)  ; 
blk=zeros(num,2) ; 

for  ii=l:num; 

OM=order(ii,l) ;0h=order(ii,2) ; 
nM=0M*(0h+l) ;nh=0h; 
blk(ii, :)=[nM  nh] ; 
name=names(ii , 1 :namelen(ii) ) ; 

eval( [' c=kron(Mmat (1 : OM+1 , 1 :0M+1) ,hmat (1 : Oh+1 , 1 : Oh+1) )*P0LY_ ' , , . . 
name  ' ; ' ] ) 

eval(['LFT_'  name  '=pss2sys( [zeros(nM,nh+l)  eye(nM) ; ' , . . . 

' zeros (nh,nM+l)  eye(nh) ; c(nM+nh+l :-l :!)''] ,nM+nh) ; ’] ) 
eval ( [ '  [LFT_ ’  name  '  ,nrem]=lftmin(LFT_'  ncime  ' ,  [nM  nh]  ,  1 , 1 , 1  ,tol)  ; ']  ) 
blk(ii, :)=blk(ii, :)-[nrem  0]; 
end 

VItVItVI'ltIt/X  LFT  construction  for  1ft ab 

LFT_ l=sbs (LFT.Xu , LFT.Xa , LFT.Xq , LFT.Xt , LFT.Xde) ; 

LFT_2=sbs (LFT.Zu , LFT.Za , LFT.Zq , LFT.Zt , LFT.Zde) ; 

LFT_3=sbs (LFT.Mu , LFT.Ma , LFT.Mq , LFT.Mt , LFT.Mde) ; 

LFT_4=sbs (0,0, 1,0,0) ; 

LFT_AB=abv(LFT_l ,LFT_2,LFT_3,LFT_4) ; 

y,  Separate  the  M  and  h  "states" 
he=cumsum(sum(blkO  0  ;me=he-blk(:  ,2)  ; 
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hs=he-blk(: ,2)+ones(num,l) ;ms=me-blk( : ,l)+ones(num,l) ; 

mindx=  []  ; 
hin<ix=  []  ; 
for  i=l:niun, 

mindx= [mindx  ms(i):me(i)]; 
hindx=[hindx  hs(i):he(i)]; 
end 

indx= [mindx  hindx] ; 

LFT_AB=reordsys(LFT_AB, indx) ; 

'/,  Reduce  the  delta  blocks  one  at  a  time 
tblk=sum(blk) 

[LFT_AB,nrem]  =  lftmin(LFT_AB,tblk,l,5,4,tol) ; 
tblk  =  tblk  -  [nrem  0] 

[LFT_AB , nrem]  =  lftmin(LFT_AB,tblk,2,5,4,tol) ; 
tblk  =  tblk  -  [0  nrem] 

'/,  Now  try  a  schur  model  reduction 

[LFT_AB,nrem]  =  If tmin2(LFT_AB, tblk, 1,5,4, tol) ; 
tblk  =  tblk  -  [nrem  0] 

[LFT_AB , nrem]  =  lftmin2(LFT_AB,tblk,2,5,4,tol) ; 
tblk  =  tblk  -  [0  nrem] 


minf o(LFT_AB) 

Clean  up  the  workspace 

clear  K  c  Mmat  hmat  ms  hs  me  he  indx  kvec  DM  Oh  blk 
clear  LFT_1  LFT_2  LFT_3  LFT_4  mindx  hindx  flag 
clear  nM  nh  neune  neimes  namelen  maxlen  ii  num  nrem 

y,  end  of  f  18poly21ft  .m 
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C.5  f  18sppoly21ft  .m;  Convert  Short  Period  Polynomials  to  LFT 


0/ 0/1/1/ 0/ 1/1/1/ t/o/vi/ o/«/ 1/1/ o/i/vo/o/./Vo/o/vvVVVVVVVVVVVVVVVVVVVVVVVVW 


fl8sppoly21ft.m 

This  script  file  defines  the  parameter  normalization  factors, 
calculates  polynomial  coefficients  in  the  normalized  deltas, 
eind  constructs  LFTs  (in  deltas)  for  Za,  Zde,  Ma,  Mq,  Mde; 
these  LFTs  are  combined  to  form  lftab=[Za  1  Zde;Ma  Mq  Mde] 

Inputs:  names  =  a  string  containing  the  names  of  the  plant  parameters 
(Za,  Zq,  etc) 

POLY.xx  =  poly,  coef.s  for  the  parameters  given  in  names,  i.e. 

POLY.Za,  POLY.Zq,  etc. 

maxfit  =  max.  order  of  the  poly's  defined  in  P0LY_xx 
order  =  matrix  defining  the  order  (in  M  and  h)  of  each  poly,  for 
each  of  the  derivatives  (i.e.  Za,  Zde,  ...) 

Each  row  of  order  coeresponds  to  one  of  the  derivatives, 
eind  the  are  listed  in  the  order  that  they  are  given  in 
nzimes 

Outputs:  LFTs  for  each  of  the  derivatives  aind  the  A  and  B  matrices. 


Written  by:  Capt  Mark  Spillman  aind  Lawton  Lee 
WL/FIGC-3  Bldg  146,  WPAFB 


'/,  Modified  by:  Capt  Martin  Breton,  AFIT  (generalized  +  schur  redux) 


•/,  DEFINE  inputs 

names='Za  Zq  Zde  Ma  Mq  Mde'; 

f lag=order(2 , 1) ; 
if  flag  ==  0, 

order= [order (1 , :) ;order(3:6, :)]  ; 
names='Za  Zde  Ma  Mq  Mde'; 

LFT_Zq=POLY_Zq; 

end 

tol=[]  ; 

'/X'ltl'ItltVLVh  Parameter  normalization  factors  and  shortcut  matrices  to  convert 
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polynomial  coefficients  in  parm  to  coefficients  in  delta_parm 

'/,  Design  and  spread  Machs 
Mbar= (Mmax+Mmin) /2 ; 

Mtil=Mmax-Mbar ; 

'/,  Design  and  spread  altitudes 
hbar  =  500*(hmax+hmin) ; 
htil  =  1000*hm2LX-hbar; 

'/,  Create  the  shortcut  matrices 
Mmat=zeros(meOcfit+l)  ; 
hmat=zeros(maLxfit+l)  ; 

K=zeros(meixf  it+1)  ;kvec=l  ;K(1 ,1)=1 ; 
for  ii=l :maxf it+1 

Mmat=Mmat  +  diag((Mtil,“(0:mzixfit-ii+l))*(Mbar“(ii-l))  ,ii-l)  ; 
hmat=hmat  +  diag((htil."(0:maxfit-ii+l))*(hbar"(ii-l)) ,ii-l) ; 
if  ii<=maxfit 
kvec=conv(kvec, [1  -1]); 

K(l:ii+l,ii+l)=abs(kvec’) ; 
end 
end 

Mmat=Mmat . *K ; 
hmat=hmat . *K ; 

VltVIXVI'ItVIt  LFT  contruction  for  Za,  Zde,  Ma,  Mq,  Mde,  qbar,  Vt 

[num, names, namelen,maLxlen]  =  namstuff  (names)  ; 
blk=zeros(num,2) ; 

for  ii=l:num; 

OM=order(ii,l) ;0h=order(ii,2) ; 
nM=0M*(0h+l);nh=0h; 
blk(ii, :)=[nM  nh] ; 
ncime=ncLmes(ii,l:neimelen(ii)) ; 

eval( [' c=kron(Mmat (1 :0M+1, 1 :0M+1) ,hmat (1 :0h+l , 1 :0h+l))*P0LY_' , . . . 
name  ' ; ' ] ) 

eval(['LFT_'  name  '=pss2sys( [zeros(nM,nh+l)  eye(nM) ; ' , . , . 

’ zeros (nh,nM+l)  eye(nh) ; c(nM+nh+l :-l : 1) ' ,nM+nh) ; '] ) 
eval ( [ ' [LFT_ ’  name  ' ,nrem]=lftmin(LFT_'  name  [nM  nh] ,1,1,1, tol) ; ']) 
blk(ii, :)=blk(ii, :)-[nrem  0]; 
end 
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'IXVl'ItVItVIX  LFT  construction  for  1ft ab 

LFT.t  op=sbs (LFT.Za , LFT.Zq , LFT.Zde) ; 

LFT_bot =sbs (LFT.Ma , LFT.Mq , LFT.Mde) ; 

LFT_AB=abv (LFT.t op, LFT.bot) ; 

'/,  Separate  the  M  and  h  "states" 

he=cumsum(sum(blkO  ')  ;me=he-blk(:  ,2)  ; 

hs=he-blk(:  ,2)+ones(nuin,l)  ;ms=me-blk(:  ,l)+ones(num,l)  ; 

mindx=  []  ; 
hindx=[]  ; 
for  i=l:num, 

mindx=[mindx  ms(i):nie(i)]; 
hindx=Chindx  hs(i) :he(i)] ; 
end 

indx= [mindx  hindx] ; 

LFT_AB=reordsys(LFT_AB, indx) ; 

'/,  Reduce  the  delta  blocks  one  at  a  time 
tblk=sum(blk) 

[LFT_AB,nrem]  =  lftmin(LFT_AB,tblk, 1 ,3,2,tol) ; 
tblk  =  tblk  -  [nrem  0] ; 

[LFT.AB , nrem]  =  lftmin(LFT_AB,tblk,2,3,2,tol) ; 
tblk  =  tblk  -  [0  nrem] 

'/,  Now  try  a  schur  model  reduction 

[LFT.AB , nrem]  =  If tmin2 (LFT.AB, tblk, 1,3,2, tol) ; 
tblk  =  tblk  -  [nrem  0] 

[LFT.AB , nrem]  =  If tmin2 (LFT.AB, tblk, 2, 3, 2, tol) ; 
tblk  =  tblk  -  [0  nrem] 

minfo (LFT.AB) 

VIXXVIXVIXL  Clean  up  the  workspace 

clear  K  c  Mmat  hmat  ms  hs  me  he  indx  kvec  OM  Oh  blk 
clear  LFT.top  LFT.bot 
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clear  nM  nh  neime  names  nzimelen  maxlen  ii  num  nrem 


'/,  end  of  f  18sppoly21ft  .m 


•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/«/•/•/ •/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/»/»/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•) 


C.5.1  If  train,  m;  Minimize  the  LFT  Realization.  Note:  this  is  a  subpro¬ 
gram  called  by  both  f  18poly21ft  .m  and  f  18sppoly21ft  .m. 


•/•/•/•/•/•/•/•/•/•/•/«/ •/•/•/•/VVVVVVVVVVVVVVVVVVVVVVV«/VVVVVVVVVV*j 


'/,  If  train,  m 

function  [LFTnew,nrera]  =  lftrain(LFTbig,blk,pnum,ninputs,noutputs,tol) ; 

'/,  function  [LFTnew,nrera]  =  lftmin(LFTbig,blk,pnura,ninputs,noutputs,tol)  ; 

•/. 

'/,  Minimize  the  block  structure  of  an  LFT  by  truncating 

'/,  "uncontrollable"  &  "unobvservable"  scalar  deltas 

y,  Inputs:  LFTbig  The  LFT  to  reduce 

y,  blk  the  original  block  structure 

y,  pnum  the  choice  of  parameter  to  be  reduced 

y,  ninputs  number  of  external  inputs,  default  =  1 

y,  noutputs  number  of  external  outputs,  default  =  1 

y,  tol  tolereince  for  minimal2.m 

y,  Outputs:  LFTnew  reduced  LFT 

y,  nrem  number  of  removed  states 

if  nargin  ==  3 
ninputs  =  1; 
noutputs  =  1; 
end 

n  =  sum(blk) ; 
nd  =  blk (pnum) ; 

[a,b,c,d]  =  minf 0 (LFTbig) ; 
if  n  ~=  d 

error( 'block  structure  does  not  add  up'); 
end 
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if  (noutputs  ~=  b  |  ninputs  "=  c) 
errorC'the  input-output  structure  does  not  add  up'); 
end 

if  a  ~=  'syst' 

errorC'lft  must  be  a  system  matrix'); 
end 

y,  Re-order  the  delta_i's  so  that  delta.pnum  comes  first 
if  pnum  >  1 

nold  =  sum(blk(l :pnum-l)) ; 

LFTmid  =  reordsys(LFTbig, [nold+1 :n,  l:nold]); 
else 

LFTmid  =  LFTbig; 
end 

y,  Partition  &  unpack  the  LFT  realization,  isolating  delta.pnum 

pss  =  sys2pss (LFTmid) ; 

a  =  pss(l:nd,l:nd) ; 

b  =  pss(l:nd,nd+l:n+ninputs) ; 

c  =  pss(nd+l:n+noutputs,l:nd) ; 

d  =  pss(nd+l :n+noutputs,nd+l :n+ninputs) ; 

y,  Remove  "uncontrollable,  unobservable"  delta  channels 
if  isempty(tol)==0 

[2un,bm,cm,dm]  =  minimal2(a,b,c,d,tol)  ; 
else 

[cun,bm,cm,dm]  =  minimal2(a,b,c,d) ; 
end 

ndnew  =  max (size (am) ) ; 
nrem  =  nd  -  ndnew; 

disp( [int2str(nrem)  '  state(s)  removed']); 

y,ashift=  am+eye  (ndnew)  ;  '/,  shift  to  handle  poles  at  origin 
y,[aim,bm,cm,dm,errbnd,hsv]  =  schmr(ashift,bm,cm,dm,2,le-3) ; 
y.ndnewer  =  max(size(am)) ; 
y,am=  am  -  eye(ndnewer) ;  '/,  shift  back 

y,nrem=  nd  -  ndnewer; 

if  nrem  >  0 
pss  =  [am,bm;cm,dm] ; 
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LFTnew  =  pss2sys(pss,m2Lx(size(pss))-max(ninputs,noutputs)) ; 

y.  Shuffle  the  deltas  back  to  their  original  order 
if  pnum  >  1 

LFTnew  =  reordsys(LFTnew, [n-nold-nrem+1 :n-nrem,  1 :n-nold-nrem] ) ; 
end 
else 

LFTnew  =  LFTbig; 
end 


y,  end  of  Iftmin.m 


C.6  flSol.m;  Convert  Full  Longitudinal  LET  to  Design  Model 


•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/VVVVVVVVVVVVVVVVVVVVVVVV'!/ 


y.  flSol.m 

function  [P , Pn , plant , qmax , err , Wp , Wh , Wr , Wu , Wn , Wc] =f 18ol (LFT_AB , tblk , qmax , err , Wp , Wu , V 


y.  fl8ol 

y.  [P ,  Pn ,  plant ,  qmax ,  err ,  Wh ,  Wp ,  Wr ,  Wu ,  Wn]  =  f  18ol  (LFT.AB ,  tblk ,  qmax ,  err) 

y. 


y,  This  function  forms  the  open-loop  interconnection  for  an  F-18 
y,  full  longitudinal  H_inf  LPV  design  problem,  using  LFT  parameter  dependence, 
y,  The  user  can  specify  the  maximum  pitch  rate,  maximum  tracking  error,  and 
y,  commeuid,  ideal  and  performcince  models  if  desired. 

y. 


y,  Inputs:  LFT.AB  Reduced  LFT  for  [A,B]  (system  matrix) 

y,  tblk  [size  of  parameterl  size  of  parcimeter2] 

y,  qmeix  Maximum  pitch  rate  (deg/sec) 

y,  err  Maximum  pitch  rate  tracking  error 

y,  example:  if  err=10,  majc  error  <  .l*qm2LX 

y. 


y,  Outputs: 

y. 

y. 

y. 


P  weighted  open-loop  interconnection  w/  uncertainty  (sys  matrix) 

Pn  weighted  open-loop  ic  w/o  uncertainty  (sys  matrix) 

NOTE:  both  P  eind  Pn  are  LPV 
plaint  Plant  only  (sys  matrix) 
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'/,  qmax  Maximum  pitch  rate  (deg/sec) 

'/,  err  Maximum  pitch  rate  tracking  error  (’/.qrnax) 

'/,  W's  All  the  weightings 

•/. 

if  nargin  ==  0 

dispC'Usage;  [P,Pn,plcint, qmax, err, Wp,Wh,Wr,Wu,Wn,Wc]  = 

'f 18ol(LFT_AB,tblk, qmax, err, Wp,Wu,Wc) ; ’ ) 
return; 
end 

if  nargin<7 ; Wc= [] ; end 
if  nargin<6 ; Wu= [] ; end 
if  nargin<5 ; Wp= [] ; end 


'/,  Command  Weight 

y,  Need  to  be  careful  putting  dynamic  filters  here.  They  tend  to  cause  a 

y,  an  unacceptable  time  delay  in  the  pitch  time  response. 

y,wnl=4; 

y,Wr=nd2sys(wnl"2,conv([l  wnl]  ,  [1  wnl] )  ,qmax*pi/180)  ; 
y,Wr=qmax*pi/ 180 ; 

Wr=l; 

y,  Performance  weight 
if  Wp==[] 

Wp=nd2sys(Cl  50]  ,  [1  4]  ,4/50/(err/100*qmax*pi/180)) Less  than  err'/,  ss  error 
%Wp=nd2sys( [1  50] ,  [1  4] ,4/50* . 1745) ; 
yWp=err/100*qmax*pi/180; 
end 

•/.  Ideal  Model 

'/,  Overshoot  and  ss  error  are  less  critical  than  the  speed  of  the  response. 

zeta=.5; 

wn2=4; 

Wh=nd2sys(wn2“2, [1  2*zeta*wn2  wn2"2]); 

*/,  Noise  weight  Wn(s) 

'/,  The  small  amount  of  noise  below  was  added  to  make  the  weighted  plant  regular 

sig_a  =  .01;  '/,  AOA  noise  variance  (rad),  realistic  value  =.01 

sig_q  =  .01;  '/,  rate  gyro  variance  (rad),  realistic  value  =.03 

Wn  =  daug(sig_a,sig_q) ; 
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'/,  Actuator  model 

au=20.2;  '/.  Elevator  bandwidth 

Act=pck(-au,au,  [l;-au] ,  [0;au]) ;  '/,  output  =  [u;du/dt] 


y,  Control  Weight 
if  Wc==[] 

Wc=l/70; 

end 

y,  LPV  plcint  with  LFT  parameter  dependence 
[ns , nx , nxu , nO]  =  minf o (LFT_AB) ; 
nu  =  nxu  -  nx; 

[dOO,blft,bO,dlft]  =  unpck(LFT_AB) ;  '/,  Unpack  LFT  realization 

aO  =  dlft( : , 1 :nx) ; 

b2  =  dlft( : ,nx+l :nxu) ; 

cO  =  blft( : , 1 :nx) ; 

d02  =  blft( : ,nx+l :nxu) ; 

y,c2  =  [0  1  0  0;0  0  1  0]  ;  y,  form  c2  and  d22 
c2=[0  1  0  0;0  0  10]; 

[ny  dum]  =  size(c2); 
d20  =  zeros(ny,n0) ; 
d22  =  zerosCny ,nu) ; 

plcint  =  pck(a0,[b0  b2]  ,  [c0;c2]  ,  [dOO  d02;d20  d22])  ; 

y,  Dummy  is  used  to  make  sure  the  number  of  controlled  outputs  equals 
y,  the  number  of  exogenous  inputs 
dummy=0 ; 


msize=tblk(l) ; 
hsize=tblk(2) ; 
delsize=msize+hsize; 
aoa=delsize+l; 
q=delsize+2; 

y,  Weighted  plant  without  uncertainty 

systemnames  =  'plcint  Act  Wp  Wc  Wh  Wr  Wn  dummy'; 

inputvar  =  ' [delm(5) ;delh(5) ;ref_cmd;noise(2) ;cont_cmd] ' ; 

outputvar  =  ' [plant(l:10) ;Wp;Wc;dummy;Wr;plant(ll:12)+Wn] ' ; 

input _to_Act  =  ' [cont.cmd] ' ; 

input _to_Wr  =  ' [ref _cmd] ' ; 
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input_to_Wc  =  '[Act(2)]'; 

input_to_Wp  =  ' [Wh-plant(12)] ' ; 

input _to_plant  =  ’ [delm;delh;Act(l)] ’ ; 

input _to_Wh  =  ' [Wr] ’ ; 

input_to_Wn  =  ’ [noise] ' ; 

input _to_duinmy  =  '  [ref_cmd]  ' ; 

cleanupsysic  =  'yes'; 

sysoutnajne  =  'fl8ic_n' ; 

sysic; 

Pn=f 18ic_n; 

y.  Uncertainty  weight 
if  Wu==[] 

Wu  =  nd2sys([l  40] ,  [1  10000] , 10000/40*le-2) ; 
end 


'/,  Weighted  plant  with  input  additive  uncertainty 

systemnames  =  'plant  Act  Wp  Wc  Wh  Wr  Wu  Wn  dummy'; 

inputvar  =  ' [delm(5) ;delh(5) ;unc_in;ref_cmd;noise(2) ;cont_cmd] ' ; 

outputvar  =  ' [plant (1 : 10) ;Wu;Wp;Wc;dummy;Wr;plant (11 : 12)+Wn] ' ; 

input_to_Wr  =  ' [ref_cmd] ' ; 

input _to_Act  =  ' [cont_cmd] ' ; 

input _to_Wc  =  ' [Act (2)]  ' ; 

input_to_Wp  =  ' [Wh-pl2mt(12)] ' ; 

input _to .plant  =  ' [delm;delh;Act(l)+unc_in] ' ; 

input_to_Wh  =  ' [Wr] ' ; 

input_to_Wu  =  '[Act(l)]'; 

input_to_Wn  =  ' [noise] ' ; 

input _to .dummy  =  ' [ref .cmd] ' ; 

cleanupsysic  =  'yes'; 

sysoutname  =  'fl8ic' ; 

sysic; 

P=fl8ic; 


y,  end  fl8ol.m 


«/•/•/•/ •/•/•/o/o/vvo/ •/•/•/«/•/•/•/»/»/•/•/•/•/•/•/•/»/•/ •/•/»/•/•/•/»/•/•/•/•/•/•/•/•/•/•/•/“/•/•/«/•> 


C.7  fl8olsp.m;  Convert  Short  Period  LFT  to  Design  Model 
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'/,  flSolsp.m 

f unct  ion  [P ,  Pn ,  plant ,  qmcix ,  err ,  Wp ,  Wh ,  Wr ,  Wu ,  Wn , Wc]  =f  ISolsp  (LFT_AB ,  tblk ,  qmax ,  err , Wp , Wi 


'/,  flSol 

'/,  [P ,  Pn ,  plant ,  qmax ,  err ,  Wh ,  Wp ,  Wr ,  Wu ,  Wn]  =  f  ISolsp  (LFT_AB ,  tblk ,  qmax ,  err) 

•/. 


'/,  This  function  forms  the  open-loop  interconnection  for  an  F-18 
'/,  full  longitudinal  H_inf  LPV  design  problem,  using  LFT  parameter  dependence. 
'/,  The  user  can  specify  the  meiximum  pitch  rate,  meiximum  tracking  error,  and 
'/,  commeind,  ideal  and  performance  models  if  desired. 

•/. 

’/,  Inputs:  LFT.AB  Reduced  LFT  for  [A,B]  (system  matrix) 

y,  tblk  [size  of  parameterl  size  of  parameter2] 

y,  qmax  Maximum  pitch  rate  (deg/ sec) 

y,  err  Maximum  pitch  rate  tracking  error 

y,  example:  if  err=10,  max  error  <  .l*qmax 


y. 

y,  Outputs:  P  weighted  open-loop  interconnection  w/  uncertainty  (sys  matrix) 

%  Pn  weighted  open-loop  ic  w/o  uncertainty  (sys  matrix) 

%  NOTE:  both  P  and  Pn  are  LPV 

y,  plaint  Plant  only  (sys  matrix) 

y,  qmax  Maximum  pitch  rate  (deg/sec) 

y,  err  Maximum  pitch  rate  tracking  error  ('/.qmeix) 

y,  W’s  All  the  weightings 

y. 


if  nargin  ==  0 

disp('Usage:  [P,Pn, plant, qmax, err, Wp,Wh,Wr,Wu,Wn,Wc]  =  ',••• 
'f  18olsp(LFT_AB, tblk, qmaix, err, Wp,Wu,Wc)  ; ') 
return; 
end 


if  nargin<7 ; Wc= [] ; end 
if  nargin<6;Wu=[] ;end 
if  nargin<5 ; Wp= [] ; end 


y,  Commaind  Weight 

y.  Need  to  be  careful  putting  dynamic  filters  here.  They  tend  to  cause  a 


'/,  am  unacceptable  time  delay  in  the  pitch  time  response. 
'/.wnl=4; 

'/,Wr=nd2sys(wnl~2,conv(  [1  wnl]  ,  [1  wnl]  )  ,qmax*pi/180)  ; 
'/,Wr=qmax*pi/ 180 ; 

Wr=l; 


'/,  Performance  weight 
if  Wp==[] 

Wp=nd2sys([l  50]  ,  [1  4]  ,4/50/(err/100*qmax*pi/180)) Less  than  err'/,  ss  error 
5CWp=nd2sys(  [1  50]  ,  [1  4]  ,4/50* .  1745)  ; 

'/,Wp=err/  100*qmax*pi/180 ; 
end 

'/,  Ideal  Model 

'/,  Overshoot  and  ss  error  are  less  critical  than  the  speed  of  the  response. 

zeta=.5; 

wn2=4 ; 

Wh=nd2sys(wn2~2,  [1  2*zeta*wn2  wn2“2]); 

'/,  Noise  weight  Wn(s) 

'/,  The  small  amount  of  noise  below  was  added  to  mahe  the  weighted  plant  regular 

sig_a  =  .01;  %  AOA  noise  variance  (rad),  realistic  value  =.01 

sig_q  =  .01;  V,  rate  gyro  variance  (rad),  realistic  value  =.03 

Wn  =  daug(sig_a,sig_q) ; 

'/,  Actuator  model 

au=20.2;  %  Elevator  bandwidth 

Act=pck(-au,au, [l;-au] , [0;au]) ;  %  ouput  =  [u;du/dt] 


'/,  Control  Weight 
if  Wc==[] 

Wc=l/70; 

end 

*/,  LPV  plant  with  LFT  parameter  dependence 
[ns,nx,nxu,n0]  =  minf o(LFT_AB) ; 
nu  =  nxu  -  nx; 

[dOO  ,blft  ,b0  ,dlft]  =  unpck(LFT_AB)  ;  '/,  Unpack  LFT  realization 
aO  =  dlft( : , 1 :nx) ; 
b2  =  dlft( : ,nx+l :nxu) ; 
cO  =  blft( : , 1 :nx) ; 
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d02  =  blft( : ,nx+l rnxu) ; 

c2  =  eye(2);  '/,  form  c2  and  d22 

[ny  dum]  =  size(c2) ; 

d20  =  zeros(ny,nO) ; 

d22  =  zeros(ny,nu) ; 

plant  =  pck(aO,[bO  b2] , [cO ; c2] , [dOO  d02;d20  d22] ) ; 

y,  Dummy  is  used  to  make  sure  the  number  of  controlled  outputs  equals 
y,  the  number  of  exogenous  inputs 
dummy =0 ; 

msize=tblk(l) ; 
hsize=tblk(2) ; 
delsize=msize+hsize ; 
aoa=delsize+l ; 
q=delsize+2; 

y,  Weighted  plant  without  uncertainty 

systemnajnes  =  'pleint  Act  Wp  Wc  Wh  Wr  Wn  dummy'; 

inputvar  =  ' [delm(3) ;delh(3) ;ref_cmd;noise(2) ;cont_cmd] ' ; 

outputvar  =  ' [plant (1:6) ;Wp;Wc ; dummy ;Wr; plant (7:8) +Wn] ' ; 

input_to_Act  =  ' [cont.cmd] ' ; 

input _to_Wr  =  ' [ref _cmd] ' ; 

input _to_Wc  =  ' [Act (2)] ' ; 

input _to_Wp  =  ’ [Wh-plant(8)] ' ; 

input_to_pleait  =  '  [delm;delh;Act(l)]  ' ; 

input _to_Wh  =  ' [Wr] ' ; 

input _to_Wn  =  ' [noise] ' ; 

input_to_dummy  =  ' [ref_cmd] ' ; 

cleeinupsysic  =  'yes'; 

sysoutname  =  'fl8ic_n'; 

sysic; 

Pn=f 18ic_n; 

y,  Uncertainty  weight 
if  Wu==  [] 

Wu  =  nd2sys([l  40] , [1  10000] , 10000/40*le-2) ; 
end 

y,  Weighted  plant  with  input  additive  uncertainty 

systemnames  =  'plant  Act  Wp  Wc  Wh  Wr  Wu  Wn  dummy'; 

inputvar  =  ' [delm(3) ;delh(3) ;unc_in;ref _cmd;noise(2) ;cont_cmd] ' ; 
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outputvar  =  '  [plant (1:6)  ;Wu;Wp ;Wc;duinmy;Wr ; plant (7:8) +Wn]  ' ; 

input _to_Wr  =  ' [ref _cmd] ' ; 

input _to_Act  =  ’  [cont_cnid]  ' ; 

input _t 0 _Wc  =  '[Act(2)]'; 

input _to_Wp  =  ' [Wh-plamt(8)] ' ; 

input _to_plant  =  ' [delm;delh;Act(l)+unc_in] ' ; 

input _to_Wh  =  ’ [Wr] ‘ ; 

input _t o _Wu  =  '[Act(l)]'; 

input _to_Wn  =  ’ [noise] ' ; 

input _to_duinmy  =  ’  [ref _cmd]  ’  ; 

cleanupsysic  =  ’yes’; 

sysoutname  =  ’fl8ic’ ; 

sysic; 

P=f 18ic; 


'/,  end  of  fl8olsp.m 


•/  •/  V  V  •/  •/  •/  V  V  5i5(y,y,5(y  V 


^0  /D  /•  /•  /•  n  !%  !%  !%  /•  /a  /a 


C.8  dkdM  1 8. m;  Per/orm  Iterations 


«/i/»/ 1/0/ •/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/ VVVVVVVVVVVVVVVW'/'/VVVVVVVVVVVVVVVVVVVVVVVVVV!/ 


y,  dkdM18.ni 

'/,  this  is  a  modif ication/extension  of  MATLAB’s  dkit 

•/.DEFINE  SOME  THINGS  BELOW 
clear 

format  short  e; 

fl8poly  '/,  for  full  longitudinal  design 

fl8poly21ft 

w=logspace(-6,6, 100) ; 

'/,fl8sppoly  '/,  for  short  period  design 
'/,f  18sppoly21ft 

*/,  Pick  some  weights 
Wp=nd2sys([l  100] , [1  4],. 4/10. 2); 

'/.Wp=nd2sys([l  10000]  ,  [1  4]  ,4/10000*10)  ; 
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'/.Wp=10; 

'/,Wp=nd2sys(  [1  4]  ,  [1  50]  ,50/4*20)  ; 

Wu=nd2sys([l  100] , [1  10000] , 10000/100*le-3) ; 

•/.Wc=l; 

Wc=l/70; 

Wn=.01; 

'/,  Plot  the  frequency  responses  of  the  weights 
Wpfrsp=frsp(Wp,w) ; 

Wufrsp=frsp(Wu,w) ; 

Wcfrsp=frsp(Wc,w) ; 

Wnfrsp=frsp(Wn,w) ; 

vplot('liv,lm' jWpfrsp, ’r-' ,Wufrsp, 'w-' ,Wcfrsp, 'b-' ,Wnfrsp, 'g-0 
title( 'Weights') 

legendC'r-' , 'Wp' , 'w-' , 'Wu' , 'b-' , 'Wc' , 'g-' , 'Wn') 

xlabel ( ' Frequency  (rad/sec) ' ) 

ylabel ( ' Magnitude ' ) 

axis([  le-1  le6  le-4  le2] ) 

pause (1) 

'/,  Form  the  weighted  plant 

[P,Pn,plant ,qmctx,err ,Wp,Wh,Wr,Wu,Wn,Wc]  =  f 18ol(LFT_AB,tblk, 1 ,15,Wp,Wu,Wc) ; 
*/,CP,Pn, plant  ,qmax,err,Wp,Wh,Wr,Wu,Wn,Wc]  =  f  18olsp(LFT_AB,tblk,  1 , 15,Wp,Wu,Wc)  ; 

msize=tblk(l) ;  Y,  tblk  defined  in  fl8poly21ft  call 

hsize=tblk(2) ; 

delsize=msize+hsize ; 

aoa=delsize+l ; 

q=delsize+2; 

THE  FOLLOWING  DEFINES  THE  DIMENSIONS  OF  THE  SYSTEM  YX/JJ. 

nmtvb=4;  */,  #  of  measured  time  varying  blocks  (2  parameters  for  both  P,K) 
ntvb=nmtvb;  '/.  #  of  time  varying  blocks 

y,  NOTE:  in  general  ntvb<=nmtvb,  but  currently  let  be  = 

y,  define  dimension  of  paraoneter  blocks  (theata.p  aind  theata.k) 

nv_p=delsize; 

nz_p=nv_p ; 

nv_k=nv_p ; 

nz_k=nv_k; 
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'/,  define  dimension  of  uncertain  block  (delta) 
nv_u=l ; 
nz_u=l ; 

'/,  define  dimension  of  performance  block  (delta.perf) 

nd=3; 

ne=3; 

'/,  define  dimension  of  the  "performance"  block  used  for  Imi 

y,  this  includes  d  to  e  plus  any  unmeasured  uncertainty  (nv_u  to  nz_u) 

nw=nd+nv_u;  V,  currently  assume  ne+nz_u=nd+nv_u ; 

y,  define  dimension  of  measures  and  actuation  of  P 

ny_p=3; 

nu_p=l ; 

y,  calculate  the  number  of  I/O's  of  the  controller 
ny=ny_p+nv_k ;  */,  #  inputs  to  K 
nu=nu_p+nz_k ;  */,  #  outputs  from  K 


y,  define  block  strunctures 
blk_p=[msize,0;hsize,0] ; 
blk_k=blk_p; 
blk_pk= [blk_p ; blk_k] ; 

blk_wu=[blk_k;blk_p;nv_u,nz_u;nd,ne]  ;  */,  blk  with  uncertainty  and  d2e 

y.  NOW  GET  SOME  BLOCKS  FOR  A  REORDERED  PARAMETER  BLOCK 
y,  NOTE:  THESE  ARE  USED  FOR  IN  A  MU  CALCULATION  FOR  COMPARISON 
blk_pkro=blk_p*2 ;  */,  blk.pk  reordered  (i.e.  after  par.rord.m  is  run) 
blk_wuro=[-blk_pkro;nv_u,nz_u;nd,ne]  ;  '/,  blk  with  uncertainty  and  d2e,  used  w/mu 
y,  for  REAL  parameters 

'/i  remove  any  zero  dimension  blks  in  the  blocks 
z_ind=f ind(blk_wu( : ,1)==0); 
blk_wu(z_ind, : )= [] ; 

z_ind=find(blk_wuro( : ,1)==0) ; 
blk_wuro(z_ind, :)=[] ; 
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set  some  options 

key_opt=0;  '/,  if  1  stops  for  input  from  keyboard  once  every  iteration 

y////////////////////.y.%%y//.y.y.*/.y.ny.y.  end  of  user  inputs  y.y.y. 


y//////////////////j//.y.y////.y.y;/.y//////.y.y//.y.y.y. 


CALCULATE  SOME  STUFF  llVhl 


nv=nv_p+nv_k+nv_u ; 
nz=nz_p+nz_k+nz_u ; 


nim=nv+nd;  '/,  #  in  M,  where  M=closed  loop 
nom=nz+ne;  '/,  #  out  M 


blk_dim=dim_blk(blk_wu) ; 


nb=ynum(blk_wu) ; 

blk_c=[colsum(blk_dim(l:nmtvb, :)) ;  '/.  this  blk  is  used  w/  mu  when  const.  D’s 
blk_wu(nmtvb+l:nb, :)]  ;  '/,  are  already  wrapped  in 

y,  initilize  D’s 
Dl=eye(nom) ; 

Dr=eye(nim) ; 

Dri=eye(nim) ; 

y,  need  a  P  that  has  proper  inputs  and  outputs  so  L,  J,  and  K 
y,  can  be  folded  into  P 

y,  i.e.  weint  to  pass  parameters  from  DELTA  thru  P  to  K 
Pk=param_k(P,nz_p,nv_p) ; 

y,  now  reorder  the  parameter  channels  so  each  parameter  is  grouped 
y,  see  ">>help  par.rord"  for  more  info 
Pkro=par_rord(Pk,blk_p) ; 


if (key_opt==l) ,  keyboard;  end 


y//////.y////////////.y.y.y.y.y.y.y.y.y.y.y.y.y.y. 


START  LOOP  y.y.y.y. 


k_count=0;  ’/,  #  of  controllers  found 
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dodkit=l ; 


while (dodkit==l) ; 
if  (k_coiint>0) 

[Dsl , Dsr] =dbf it (DcMDci , blk_wu ,nmtvb , w , mmb.DcMDci , rowd.DcMDci , sns.DcMDci) ; 

'/,  similar  to  musynfit 

Dsri=minv(Dsr) ; 

Dl=iiimult  (Dsl , Del)  ; 

Dri=mmult(Dcri,Dsri) ; 

DMDi=mmult(Dl,M,Dri) ; 

DMDig=frsp(DMDi,w) ; 
msv_DMDi=sel (vsvd(DMDig) ,1,1); 

'/,  determine  what  type  of  plot  to  use 

if (pkvnorm(mmb_DcMroDci)<3  &  pkvnorm(mmb_DcMDci)<3  &  pkvnorm(msv_DMDi)<3) 
plot_type='liv,m' ; 
else 

plot_type='liv,lm' ; 
end 

vplot(plot_type,mmb_DcMroDci, 'r' ,mmb_DcMDci, 'g' ,msv_DMDi, 'b') 
t it le ( ' mu_up (DcMroDci ) =red ,  mu_up (DcMDci ) =green ,  msv (DMDi) =blue ' ) 
pos(['Q  The  red  is  the  mu  upper  bound  using  dynamic  D  scales’,...) 

’®  The  green  is  an  actual  mu  upper  bound  using  constant  and’ , . . . 

’  dynaimic  D  scales  ’ , . .  . 

’Q  The  blue  is  the  msv  of  the  closed  loop  when  TFs  are  used  in’ , . . . 

’  place  of  the  @  dynamic  D  scales  of  the  green.  ’,... 

’Q  Q  NOTE:  if  the  option  "cmp.opt"  is  set  to  zero  there  is  no  red  curve  ’]) 

'/,*/,  the  following  is  a  check  to  see  what  the  hinf  norm  of  the  weighted 
'/•%  (by  Ds,  L,  and  J)  closed  loop  system  is  using  L,  J,  and  K  from  the 
'/X  previous  step,  the  hinf  norm  after  finding  new  L,  J,  and  K  should 
'/,%  be  less  them  or  equal  to  this  norm;  worst  case  it  would  be  equal  if 
'/,*/,  the  minimization  returned  the  same  L,  J,  and  K  as  the  previous  step 
hinf norm (DMDi) 

disp(’The  hinf  norm  after  the  next  step  should  be  <=  that  shown  above’) 
'/,pos(’The  hinf  norm  after  the  next  step  should  be  <=  that  shown  above’); 
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else 

Dsl=eye(noin)  ; 

Dsri=eye(nim) ; 
end 

*/,  NOW  need  to  scale  P  so  can  iterate 
*/,  First  scale  Pk,  then  select  P 
Dslp=daug(Dsl,eye(ny)) ; 

Dsrip=daug(Dsri,eye(nu)) ; 

DsPkDsi=inmult(Dslp,Pk,Dsrip)  ; 

DsPDsi=sel(DsPkDsi ,nz_k+l :nz+ne+ny_p,nv_k+l:nv+nd+nu_p) ; 

'/•'/t*/.'/.  temp 

if  k_count==0 

[k,DwMDwi,gf]=hinfsyn(DsPDsi,ny,nu,  .2,200, le-2) ;  */,  checks  LTI  controller 
[k,DwMDwi,gf]=hinfsyn(Pk,ny,nu,  .2,1000,le-2) ;  */,  checks  unsealed  LPV 

end 

vmm,  end  temp 

'/,  CALL  LMI  FUNCTION  THAT  RETURNS  K  AND  Dc 

do  LMI  iterations  to  determine  gamma  and  find  L,J 

gamma=input(' ENTER  VALUE  OF  TARGET  GAMMA,  OR  999  FOR  GEVP  0; 
while  looptst2(geaiima,l) 

gamma=input ( 'try  again  -  ENTER  VALUE  OF  TARGET  GAMMA  OR  999  '); 

end 

if  gcimma==999 

opt=[le-2  200  le7  5  0]; 

[Ps , R , S , L , J , gamma] =nshinf lmi2 (DsPDsi , abs (blk_p) , nw , opt , 0 ) ; 
else 

dispC' BEGINNING  OF  ITERATION  LOOP'); 

gtemp=-l; 

while  gamma~=999 

disp( ['target  gamma  set  at  '  num2str(gamma)] ) ; 
opt=[0  200  le7  30  0]; 

[Ps,R,S,L,J,t] =nshinf Imi (DsPDsi , abs (blk_p) , nw , gamma , opt , - le-8) ; 
if  t<0 

gtemp=gamma; 

Ltemp=L; 

Jtemp=J ; 


0; 


Ptemp=Ps ; 
end 

gaiiima=input(' ENTER  NEW  TARGET  GAMMA,  OR  999  TO  MOVE  ON 
while  looptst2(gainma,l) 

gainma=input ( 'try  again  -  ENTER  VALUE  OF  TARGET  GAMMA  OR  999  ') 

end 
end 

if  gtemp>-l 

disp(['the  solution  for  gamma  =  '  num2str(gtemp)  '  will  be  used']); 
gcLmma=gtemp ; 

L=Ltemp; 

J=Jtemp; 

Ps=Ptemp; 

else 

error ('Problem  is  not  feasible.  Aborting  dkdit.m'); 
end 
end 


[K] =hinf  syn (Ps , nv_p+ny_p , nv_p+nu_p , 0 , gamma , 1 e-2 ) ; 


k_count=k_ count +1 ; 

eval( ['K_' ,num2str(k_count) , '=K; ']) ; 


y,  reassign  L  to  Del  aind  J  to  Deri  (D  constant  left  aind  right  inverse) 
Dcl=sqrtm(L) ; 

Dcri=sqrtm(J) ; 
y,  Dcri=Dcl\eye(2*nv_p+nw) 
y.  Dcl=chol(L); 
y,  Dcri=Dcl\eye(2*nv_p+nw) ; 

y  form  closed  loop 
M=starp(Pk,K) ; 

Mg=frsp(M,w) ; 

DcMDci*mmult (Del, M, Deri) ; 

DcMDcig=mmult (Del ,Mg,Dcri)  ;  */,  scale  closed  loop  freq  resp 

y.y.  CHECK  THE  HINF  NORM  WITH  THAT  OBTAINED  BY  LMI 

yy,  The  following  checks  the  hinf  norm  of  the  weighted  (by  Ds,  L,  and  J) 
yy  closed  loop  system  using  the  L,  J,  and  K  just  obtained  eind  the  Ds 
yy  from  the  previous  step,  this  should  be  similar  to  the  gamma  obtained 
yy  by  the  LMI  feasability  search  AND  should  be  <=  the  norm  checked  using 
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'/,%  the  L,  J,  eoid  K  from  previous  step,  which  was  checked  above 
DcDsMDsiDci=inmult(Dcl,starp(DsPkDsi,K)  ,Dcri) ; 
hinfnorm(DcDsMDsiDci) 

dispC'The  hinf  norm  of  the  new  closed  loop  system  is  shown  above’) 
'/,pos([’  The  hinf  norm  of  the  new  closed  loop  system  is  shown  above’,... 
'/,  ’0  Check  this  value  with  that  of  gfin  given  by  the  LMI’]); 

[bnd_DcMDci , rowd_DcMDci , sns.DcMDci ,rp] =mu2 (DcMDcig ,blk_c) ; 
mmb_DcMDci=sel(bnd_DcMDci,l,l) ; 
pk_mmb_DcMDci=pkvnorm(mmb_DcMDci) ; 

disp(’  ’) 

dispC’The  peak  of  the  mu  upper  bound  (green  curve)  is’); 
disp(pk_mmb_DcMDci) ; 


dodkit=input (’Press  return  to  do  real  mu,  or  1  to  skip  it.  ’); 
if  (isempty(dodkit) ) 

VI,  FORM  CLOSED  LOOP  WITH  REORDERED  SYSTEM  SO  CAN  DO  A  MU  CALC  FOR  COMPARISON 
DcMroDci=par_rord(DcMDci ,blk_p) ; 

DcMroDcig=frsp(DcMroDci,w) ; 

[bndro , d_DcMroDci , sns_DcMroDci , rp] =mu(DcMroDcig , blk_wuro) ; 

mmb_DcMroDci=sel(bndro,  1 , 1)  ;  */,  mu  upper  bound  with  full  D’s  on  theta  block 

pk_mmb_DcMroDci=pkvnorm (mmb.DcMroDci ) ; 

disp(’  ’) 

disp(’The  peak  of  the  (real  theta)  mu  upper  bound  (red  curve)  is’); 
disp(pk_mmb_DcMroDci) ; 

else 

mmb_DcMroDci=mmb_DcMDci ; 
end 


'/,  determine  what  type  of  plot  to  use 
if (pkvnorm(mmb_DcMroDci)<3  &  pkvnorm(mmb_DcMDci)<3) 
plot_type=’liv,m’ ; 
else 

plot_type=’liv,lm’ ; 
end 

vplot (plot_type,mmb_DcMroDci , ’r’ ,mmb_DcMDci, ’g’) 
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eval( ['plot_type_’ ,num2str(k_count) , '=plot_typG; '] ) ; 
eval( ['mmb.DcMroDci.' ,num2str(k_count) , '=nmib_DcMroDci; '] ) ; 
eval(  ['mmb_DcMDci_ '  ,niim2str(k_count)  ,  '=mmb_DcMDci;  '])  ; 

titleC'mmb.DcMroDci  (red),  mmb.DcMDci  (green),  NOTE:  green  is  a  true  upper  bound' 
disp(['  Q  The  green  plot  is  a  true  upper  bound  on  mu  for  this  system.',... 

'Q  The  red  plot  is  the  upper  bound  using  dynamic  scales  for  all  blocks',... 

'0  0  NOTE:  red  curve  scales  can  only  be  kept  in  last  iteration  ']) 


y//.*/.y.y.y.y.y.y.y.'/.y.y.'/.r/.y.*/.y.  do  some  analysis  y.y.y.y.y.y.y.y. 

do_m=input(' Would  you  like  to  do  an  M  analysis?  y=yes,  return=no  ','s'); 
if (~isempty(do_m)) 

'/#'/•'/•'///•  look  at  some  M  analysis  to  see  what  is  driving  the  problem 
[dl , dr] =unwrapd2 (rowd.DcMDci ,blk_c) ; 
dmdig=mmult (dl , DcMDcig , vinv (dr) ) ; 
msv_DcMDci=vnorm(dmdig) ; 
blk_m= [nv , nz ; nd , ne] ; 
msv=blknorm(dmdig,blk_m) ; 
msv_zv=sel(msv, 1 , 1) 
msv_zd=sel (msv ,1,2) 
msv_ev=sel (msv ,2,1) 
msv_ed=sel(msv,2 ,2) 


subplot (2,2, 1) ,  vplot (plot .type, mmb_DcMDci, 'r' ,msv_zv, 'g') 
title('Robust  Stability,  mu(Mzv)') 

subplot (2, 2, 2) ,  vplot (plot .type, mmb.DcMDci, 'r' ,msv.zd, 'g') 
title( 'msv(DMzdDi) ' ) 

subplot (2,2,3) ,  vplot (plot . type, mmb.DcMDci, 'r' ,msv.ev, 'g') 
title( 'msv(DMevDi) ' ) 

subplot (2,2,4) ,  vplot (plot .type, mmb.DcMDci, 'r' ,msv.ed, 'g' ) 

title ( 'Nominal  Performance,  msv(Med)') 

pos('p') 

blk.m2= [nv-nv.u , nz-nz.u ; nv.u , nz.u ; nd , ne] ; 
msv2=blknorm(dmdig,blk.m2) ; 

y,  note:  p  denotes  parameters  and  u  denotes  uncertainty 

msv.pp=sel (msv2 ,1,1); 

msv.pu=sel (msv2 ,1,2); 

msv.pd=sel (msv2 ,1,3); 

msv.up=sel (msv2 ,2,1); 

msv.uu=sel (msv2 ,2,2); 

msv.ud=sel (msv2 ,2,3); 
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insv_ep=sel(msv2,3,l) ; 
msv_eu=sel(msv2,3,2) ; 
msv_ee=sel(msv2,3,3) ; 


subplot (3,3,1),  vplot (plot_type,inmb_DcMDci ,  'r '  ,iiisv_pp,  'g' ) 
title ( 'Parameter  input’) 
ylabelC' Parameter  output’) 

subplot (3,3,2),  vplot (plot_type ,mmb_DcMDci , ’r ’ ,msv_pu, ’g’ ) 
title (’Uncertainty  input’) 

subplot (3, 3, 3) ,  vplot (plot_type,mmb_DcMDci , ’r' ,msv_pd, ’g’ ) 
t  it  1  e  ( ’  Di  sturbcince/Ref erence  input ' ) 

subplot (3,3,4) ,  vplot (plot_type ,mmb_DcMDci , ’r ’ ,msv_up , ’g ’ ) 
ylabel( ’Uncertainty  output’) 

subplot (3,3,5) ,  vplot (plot_type ,mmb_DcMDci , ’ r ’ ,msv_uu, ’ g ’ ) 
subplot (3,3,6) ,  vplot (plot_type ,mmb_DcMDci , ’r ’ ,msv_ud, ’g’ ) 
subplot (3,3,7) ,  vplot (plot_type ,mmb_DcMDci , ’r ’ ,msv_ep , ’g’ ) 
ylabel ( ’ Error/Performance  output ’ ) 

subplot (3,3,8) ,  vplot (plot_type,mmb_DcMDci , ’r ’ ,msv_eu, ’g’ ) 
subplot (3,3,9) ,  vplot (plot_type,mmb_DcMDci , ’r ’ ,msv_ed, ’g’ ) 
pos(’p’) 
subplot (1,1,1) 

end  M  analysis 

Check  convergence  '/*/,%'/*/, 
eval([’L_’ ,num2str(k_ count) , ’=L; ’]) ; 
rd=rowd_DcMDci ; 

eval([’rd_’ ,num2str(k_count) , ’=rd; ’]) ; 
if (k_count>l) 

disp ( ’ CHECKING  CONVERGENCE ’ ) ; 

eval( [’L_dif=L_’ ,num2str(k_ count) , ’-L_’ ,num2str(k_count-l) ,’;’]); 
max_L_dif=max(max(L_dif ) ) ; 

disp([’The  max  difference  between  L(i)  and  L(i-l)  is  ’ ,num2str(max_L_dif )]) ; 

[type,nr,rc,npts]=minfo(rd) ; 
eval( [’rdi=rd_’ ,num2str(k_count) ,’;’]); 
eval([’rdiml=rd_’ ,num2str(k_count-l) ,’;’]); 
for  k=l:rc; 

vplot(’liv,lm’,sel(rdi,l,k),’r’,sel(rdiml,l,k),’g’) 

title([’Dscale  #’ ,num2str(k) , ’  from  this  iteration  and  the  previous  one’]) 
pos('p’) 


end 

end 


END  OF  ANALYSIS  VI,ThVI, 


if (key_opt==l) ,  keyboard;  end 

dodkit=input (' Press  return  to  do  another  iteration,  or  1  to  quit.  '); 
if  (isempty(dodkit)) 
dodkit=l ; 
else 

dodkit=0 ; 
end 


end  't,  while 
y.  end  dkdMlS.m 


0/ 1/1/1/ 0/ 1/1/1/ o/i/t/vi/vvo/t/i/o/vo/i/vo/o/vvVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVyVVVVVVi 


C.8.1  param_k.m;  Add  Controller  Parameter  Block  to  Design  Model. 


i/i/i/i/i/i/i/i/i/i/i/i/i/i/i/i/ 1/1/1/ 1/1/1/ 1/1/1/ 1/1/1/ i/t/i/i/i/i/t/i/ 1/1/1/ 1/1/1/ 1/1/ 1/1/ 1/ 1/1/ i/i/i/ii 


y,  parajn_k.in 


y,  this  function  converts  an  LPV  plant  into  an  LPV 

y,  pleint  that  "passes"  the  same  parameters  to  the  controller 

y. 


y,  function  sys_out=param_k(sys_in,nz,nv)  ; 

y. 


y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 


where,  sys_in=LPV  plant 

where,  the  original  LPV  system  is  given  by 
LPV_orig=starp(del , sys_in) 

nz=#  of  inputs  to  the  delta  block  (i.e.  #  of  "uncertain" 
(parameter)  outputs  from  the  original  plant 
nv=#  of  outputs  from  the  delta  block  (i.e.  #  of  "uncertain" 
(parameter)  inputs  to  the  original  plant 


y,  the  new  LPV  system  is  given  by  LPV=starp(del_n,sys_out) 
y,  where,  del_n  =  [delk  I 
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•/. 

•/. 


del] 


function  sys_out=param_k(sys_in,nz,nv) ; 
Ca,b,c,d] =unpck(sys_in) ; 

[mtype , npo , npi , nx] =minf o ( sy s_in) ; 
a_new=a ; 

b_new= [zeros (nx,nv) ,b,zeros(nx,nz)] ; 

c_new= [zeros (nz,nx) ; c; zeros (nv,nx)] ; 

d_new= [zeros (nz,nv) , zeros (nz, npi) ,eye(nz) ; 
zeros(npo,nv) ,d, zeros (npo, nz) ; 
eye(nv) , zeros (nv, npi) ,zeros(nv,nz)] ; 

sys_out=pck(a_new ,b_new , c_new , d_new) ; 

y,  end  of  param_k.m 


•/ 0/ 0/ 0/ 1/ •/•/ 0/ 0/ 0/ V  •/ 1/ 0/ g/ 0/ 0/ •/ 5/ >/ 0/ 1/ 1/ i/ •/>/ 0/ 0/ 1/ •/•/ 1/ 0/ 0/ •/•/•/ 0/ 1/ 0/ •/ i/ •/ 0/ 0/ 0/ •/•/•/ 0/ >/ 0/ (i 


C.8.2  par_rord.m;  Reorder  Parameter  Block  for  Analysis. 


a/  a/  a/  a/  a/  a/  a/  a/  a/  a/  a/  a/  a/  a/  a/  a/  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  oy  oy  oy  oy  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ay  ai 


y,  par.rord.m 

y,  this  function  modifies  the  LPV  pleint  generated  in  param_k.m  which 
y,  "passes"  the  same  pareimeters  that  the  plant  varies  with  to  the  controller 

y. 

y.  The  modification  is  simply  a  reordering  of  the  "uncertainty" (/parameter)  I/O's 

y,  (i.e.  inputs /outputs)  so  that  the  final  delta  (which  includes  the  plant 

y,  and  the  controller  subdeltas)  contains  only  scalar*Identities  with  each  scalar 

y,  appearing  only  once  (i.e.  converts  del_k_p  to  del 

y,  where,  del_k_p  =  [delk  I 

y#  I  delp] 

y. 
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y,  where,  delk=diag(dl*Ikl ,d2*Ik2, .  . . )  and  delp=diag(dl>t‘Ipl,d2*Ip2, . . .) 

y,  NOTE:  assumed  here  that  delk=delp  i.e.  Ikl=Ipl  etc. 

y,  then , 

y,  del=[dl*diag(Ikl,Ipl)  I 

y,  I  d2*diag(Ik2,Ip2)  I 

y.  I  •••  ] 

y. 

y,  function  sys_out=par_rord(sys_in,nB,B_dim) 

y. 

y,  where,  sys_in=LPV  plant  generated  by  param.k.m 
y,  nB=  #  of  blocks  in  each  del  (i.e.  #  of  blocks  in  delp) 

y,  B_dim=a  vector  containing  the  dimension  of  each  block  of  each 

y,  delCe.g.  delp)  =[diml;dim2;  . .  .] 

y. 

y. 

y,  function  sys_out=pax_rord(sys_in,blk_p) 

y. 

y,  where,  sys_in=LPV  plant  generated  by  pareim.k.m 

y,  blk_p=ein  nB  by  2  matrix  that  defines  the  structure  of  the 

y,  parameter  block  of  P  (i.e.  delp) 

y. 


function  sys_out=par_rord(sys_in,nB_or_blk_p,B_dim) 
if(nargin<2  1  nargin>3) 

disp('useage:  sys_out=pax_rord(sys_in,nB,B_dim) ') ; 
return 

elseif (nargin==2) 

blk_dim=dim_blk(nB_or_blk_p) ; 
nB=ynum(nB_or_blk_p) ; 

B_dim=blk_dim( : , 1) ; 
elseif (nargin==3) 
nB=nB_or_blk_p ; 
end 

[mtype , npo , npi ,nx] =minf o (sys_in) ; 

B_sum=sum(B_dim) ; 

parm_order=[l:B_dim(l) ,B_sum+l :B_sum+B_dim(l)] ; 
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if(nB  >  1) 

for  blk_n=2:nB; 

sbdbnml=suin(B_diin(l  :blk_n-l))  ; 
sbdbn=sum(B_dim(l:blk_n)) ; 

parm_ordGr=[pann_order,sbdbnitil+l  :sbdbn,B_sum+sbdbnml+l  :B_sura+sbdbn]  ; 
end 
end 


out_order=[pann_order,B_suia*2+l  :npo]  ; 
in_order= [parm.order , B_sum*2+1 : npi]  ; 

sys_out=sel (sys_in, out_order , in_order) ; 


'/,  end  of  par.rord.m 


1/ 0/./V./ 


C.9  nshinf  Imi  .m;  Solve  LMI  Feasability  Problem 


•/•/•/•/ o/»/*/vo/#/o/o/ «/•/ 0/1/ j/i/Mo/o/vvo/o/VVVVVVVVVVVVVVVVVVVVVVVVyyVVVVVVVVy 


'/,  nshinf  Imi.  m 


f unct ion  [Ps,R,S,L,J,t] =nshinf Imi (P , blk , nw , gamma , opt .goal) 

•/. 


•/,  NSHINFLMI 

y. 

y. 

y. 


Find  the  solution  to  the  Apkarian's  scaled  Hinf  optimization 
problem.  The  problem  is  formulated  as  a  Feasibility  problem  using 
LMIs.  The  t  value  displayed  must  be  negative  in  order  for  a 
solution  to  exist  for  the  given  gamma. 


y. 

y,  usage: 

y. 

y,  inputs: 

y. 

y. 

y. 

y. 

y. 

y. 


[R,S,L,J]=nshinflmi(P,blk,nw, gamma , opt , go  al ) 

P  -  Open-loop  plant  in  packed  notation,  i.e.  P=pck(A,  ... 

[Bt  B1  B2],  [Ct;Cl;C2],  [Dtt  ...;Dlt  ...;  D2t  ...]), 
where  t  denotes  the  parameter  influence 
blk  -  nm  X  2  matrix  describing  the  structure  of  the  parameter 
block.  The  first  element  in  each  row  denotes  the  size  of 
the  subblock.  The  second  element  in  each  row  is  either 
0  (for  repeated  scalar)  or  1  (for  full) . 
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y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 


nw  -  Number  of  exogeneous  inputs  (must  equal  the  number  of 
controlled  outputs) 
gamma  -  Desired  Hinf  norm 

opt  -  Optional  vector  of  pareimeter  to  pass  to  the  LMI 

feasibility  solver.  opt=[  0  (iters)  (feas.  rad.)... 

(slow  prog,  tol)  (print)],  default=[0  200  le7  30  0].  See 
LMI  Lab  notes  for  more  details. 

goal  -  Optional  target  for  feasibility  problem,  def ault=-le-8. 

See  LMI  Lab  notes  for  more  details. 


y. 

y,  outputs:  R,  S  -  Solutions  to  the  feasibility  problem 
y,  L,  J  -  The  complete  scaling  matrices  derived  from  L3  and  J3 

y,  Ps  -  The  COMPLETE  weighted  plant  scaled  by  L  and  J 
y,  t  -  t<=0  indicates  the  target  gajnma  was  feasible 

y. 


y,  Written  by:  Capt  Mark  Spillman 
y.  HL/FIGC-3  Bldg  146,  WPAFB 

y. 


y,  Modified  by:  Capt  Martin  Breton,  AFIT 


if  nargin<l 

disp('  usage:  [Ps,R,S,L,J,t]=nshinflmi(P,blk,nw, gamma, opt, goal) ; ') ; 

disp('  0; 

return 

end 


y,  The  equations  described  below  are  teiken  from  "A  Convex  Characterization  of 
y,  Gain-Scheduled  Hinf  Controllers" 

y,  Determine  optional  inputs 
if  nargin<6;goal=-le-8;end 
if  nargin<5; opt= [0  200  le7  30  0] ;end 

y,  Get  nm  emd  change  blk  to  denote  the  structure  of  L3  and  J3 

nm=sum(blk( : ,1))  ; 

blk(:,2)=-(blk(:,2)-l); 

y,  Set  nz=nw,  a  current  restriction  of  the  code 
nz=nw ; 

y,  Unpack  P  and  spilt  apart  the  elements  in  B,  C  and  D 
[A,B,C,D]=unpck(P) ; 
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na=size (A , 1 ) ; nu=size (B , 2) -nw-nm ; ny=size (C , 1) -nz-nm ; 

Bt=B(:  ,l:nin)  ;B1=B(:  ,iiin+l  :nin+nw)  ;B2=B(:  ,nm+nw+l  :nm+nw+nu)  ; 

Ct=C(l:nm, :) ;Cl=C(nm+l :nm+nz, :) ;C2=C(nm+nz+l :nm+nz+ny, :) ; 
Dtt=D(l:nin,l:nm);Dtl=D(l:niti,  nm+ 1 :  nm+nw  );Dt2=D(l:nm,  nm+nw+ 1 :  nm+nw+nu)  ; 
D 1 t =D (nm+ 1 : nm+nz , 1 : nm) ; D 1 1=D (nm+ 1 : nm+nz , nm+ 1 : nm+nw) ; 

D12=D (nm+1 : nm+nz , nm+nw+1 : nm+nw+nu) ; 

D2t=D(nm+nz+l:nm+nz+ny,l:nm) ;D21=D(nm+nz+l :nm+nz+ny,nm+l :nm+nw) ; 
D22=D(nm+nz+l :nm+nz+ny, nm+nw+1 : nm+nw+nu) ; 

'/,  tolerances  (taken  from  Geihinet's  routines) 
macheps=mach_eps ; 
tolsing=sqrt (macheps) ; 
toleig=macheps“(2/3) ; 

y,  Make  the  appropriate  substitutions  to  use  a  part  of  Cabinet 's 
y,  routine  below 
DD12=[Dt2;D12] ;DD21=[D2t  D21] ; 

y,  The  following  was  taiken  from  Cabinet ’s  goptlmi.m: 

J^*:<C5|C3|C5|C5|C*:|C*:|c:t«*5i£5jC*:|C****5|C5|e**5|c*5|C5|C5|C********************************** 

y,  For  numerical  stability  of  the  controller  computation, 

y,  zero  the  sing,  values  of  DD12  s.t  I  I  B2  DD12“+  I  I  >  1/tolsing 

[u,s,v]=svd(DD12) ; 

abstol=max(toleig*norm(B2,l) ,tolsing*s(l , 1)) ; 
ratio=max([s;zeros(l,size(s,2))]) ./. . . 

max( [tolsing*abs(B2*v) ;abstol*ones(l,nu)]) ; 
ind2=find (ratio  <  1);  12=length(ind2) ; 

if  12  >  0,  s(: ,ind2)=zeros(nz,length(ind2)) ;  DD12=u*s*v' ;  end 
[u,s,v]=svd(DD2lO; 

abstol=max(toleig*norm(C2,l) ,tolsing*s(l , 1)) ; 
ratio=max([s;zeros(l,size(s,2))]) ./. . . 

max( [tolsing*abs(C2’*v) ;abstol*ones(l,ny)] ) ; 
ind2=find (ratio  <  1);  12=length(ind2) ; 

if  12  >  0,  s( : ,ind2)=zero s (nm+nw, length (ind2)) ;  DD21=v*s'*u’ ;  end 

y,  compute  the  outer  factors 
Nr=lnull([B2;DD12] ,0,tolsing); 
cnr=size(Nr,2) ; 

Nr=[Nr  zeros (na+nm+nz, nm+nw) ; zeros (nm+nw, cnr)  eye (nm+nw) ] ; 
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Ns=rnull([C2,DD21] ,0,tolsing) ; 
ciis=size(Ns,2)  ; 

Ns=[Ns  zeros(na+nm+nw,nm+nz) ;zeros(nm+nz,cns)  eye(nm+nz)] ; 

J^jlc  +  slcslc**********:^************************************************ 

setlmisC  []  ) ; 

'/,  Define  the  matrix  variables:  R,S,L3,J3 
R=lmivar(l, [na  1]); 

S=lmivar(l, [na  1]); 

L3=lmivar(l ,blk) ; 

[J3,nvar] =lmivar(l,blk) ; 

'/,  Display  variable  info 
disp(’  '); 

disp([ 'There  are  '  int2str(nvar)  '  standard  variables. '] ) 
dispC'  '); 


'/,  Define  the  individual  terms  in  the  first  two  equations 
ImitermC [1,0, 0,0] ,Nr) ; 

ImitermC [2, 0,0,0] ,Ns) ; 

ImitermC [1,1,1,R] ,A,1, 's') ; 

ImitermC [2, 1,1,S] ,A' ,1, 's') ; 

ImitermC [1, 2, 1,R] ,Ct,l) ; 

ImitermC  [2,2 ,1 ,S] ,Bt ' , 1) ; 

ImitermC [1,2,2, J3] , -gamma,!) ; 

ImitermC [2,2,2,L3] , -gamma,!) ; 

ImitermC [1,3,1,R] ,C1,1); 

ImitermC [2,3,1,S] ,B1' ,1) ; 

ImitermC [1,3, 3,0] , -gamma) ; 

ImitermC [2, 3, 3,0] , -gamma) ; 

ImitermC [1, 4, 1,J3] ,l,Bt') ; 

ImitermC [2,4,1 ,L3] , 1 ,Ct) ; 

ImitermC [1,4,2, J3] ,l,Dtt'); 

ImitermC  [2,4,2,L3] , 1 ,Dtt) ; 

ImitermC [1,4,3, J3] ,1 ,Dlt ' ) ; 

ImitermC [2,4,3,L3] , 1 ,Dtl) ; 

ImitermC [1,4,4, J3] , -gamma, 1) ; 

ImitermC [2,4,4,L3] , -gamma, 1) ; 

ImitermC [1,5 , 1 ,0] ,B1 ') ; 

ImitermC [2, 5, 1,0] ,C1) ; 
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lmiterm(Cl,5,2,0]  ,DtlO  ; 
lmiterm(C2,5,2,0] ,Dlt) ; 
lmiterm([l,5,3,0]  ,DllO  ; 

ImitermC [2, 5, 3,0] ,D11) ; 
lmiterm(Cl,5,5,0] , -gamma) ; 

ImitermC [2, 5, 5,0] , -gamma) ; 

y,  Define  individual  terms  in  the  last  two  equations 
lmiterm([-3,l,l,R]  ,1,1) ; 

ImitermC [-4, 1, 1,L3] ,1,1); 

ImitermC [-3, 2, 1,0]  ,1)  ; 

ImitermC [-4, 2, 1,0] ,1)  ; 

ImitermC  C-3,2,2,S]  ,1,1); 

ImitermC [-4,2,2, J3] ,1,1) ; 

lsys=getlmis; 

y,  Solve  Feasibility  problems 
[t,xfeas]=feaspClsys, opt, goal) ; 

R=dec2mat  Clsys ,xf eas ,R) ; 

S=dec2matClsys,xfeas,S) ; 

L3=dec2mat Clsys, xf eas, L3) ; 

J3=dec2mat  Clsys , xf eas , J3)  ; 

y,  Evaluate  Feasibility 
evals=evallmi Clsys, xf eas)  ; 

[Ihsl] =showlmiCevals , 1)  ; 

[lhs2]=showlmiCevals,2) ; 

[lhs3,rhs3]=showlmiCevals,3) ; 

[lhs4,rhs4]=showlmiCevals,4) ; 

y,  Display  Results 
dispC'  ') 

dispC’The  following  values  should  be  negative.') 

dispCE'The  maximum  eig  of  the  LHS  of  LMIl  is:  '  num2strCmaxCeigClhsl) ))] ) ; 

dispCE'The  maximum  eig  of  the  LHS  of  LMI2  is:  '  num2strCmaxCeigClhs2) ))] ) ; 

dispC'  ') 

dispC'The  following  values  should  be  positive.') 

dispCE'The  minimum  eig  of  the  RHS  of  LMI3  is:  '  num2strCmaxCeigCrhs3)))] ) ; 

dispCE'The  minimum  eig  of  the  RHS  of  LMI4  is:  '  num2strCmaxCeigCrhs4)))] ) ; 

dispCE'The  minimum  eig  of  R  is:  '  num2strCminCeigCR) ))] ) ; 
dispCE'The  minimum  eig  of  S  is:  '  num2strCminCeigCS)))]) ; 
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disp(['The  minimum  eig  of  L3  is:  '  num2str(min(eig(L3)))] ) ; 
disp([’The  minimum  eig  of  J3  is:  '  num2str(min(eig(J3)))] ) ; 

y,  Compute  L  and  J  (old  way) 
y,N=  (  J3\eye  (nm)  )  -L3 ; 

y,L=[-N'  N  zeros (nm,nz)  ;N'  L3  zeros(nm,nz) ;  zeros (nz,2*nm)  eye(nz)]  ; 
y,Jl=(-N'  )\eye(nm)+J3; 

y.J=[  J1  J3  zeros(nm,nz)  ;  J3  J3  zeros  (nm,  nz)  ;  zeros  (nz,2*nm)  eye(nz)]  ; 

y,  Compute  L  and  J  (new  way) 

J3p=chol(J3) ; 

J3pi=J3p\eye(nm) ; 

J3i=J3pi*J3pi’ ; 

N=(J3i-L3); 

Np=chol(-N) ; 

Npi=Np\eye(nm) ; 

Ni=Npi*Npi' ; 

Jl=J3+Ni; 

L=[-N  N  zeros (nm,nz) ;N  L3  zeros(nm,nz) ;  zeros(nz,2*nm)  eye(nz)] ; 

J=C  J1  J3  zeros(nm,nz) ;  J3  J3  zeros (nm,nz) ;zeros(nz,2*nm)  eye(nz)] ; 


y,  Form  the  elements  of  the  true  weighted  plant 
a=A;bl=[zeros(na,nm)  Bt  Bl];b2=[B2  zeros (na,nm)] ; 
cl= [zeros (nm.na) ;Ct;Cl] ; c2=[C2; zeros (nm,na)] ; 

dll=[zeros(nm,2*nm+nw) ; zeros (nm)  Dtt  Dtl ; zeros (nz,nm)  Dlt  Dll]; 
dl2= [zeros(nm,nu)  eye(nm);Dt2  zeros(nm) ;D12  zeros(nz,nm)] ; 
d21= [zeros (ny,nm)  D2t  D21;  eye(nm)  zeros (nm,nm+nw)] ; 
d22=[D22  zeros(ny,nm) ; zeros (nm,nu+nm)] ; 

y,  Scale  the  true  plant  so  that  L=J=I 

condL=cond(L) 

condJ=cond(J) 

y,F=chol(L)  ;Fi=F\eye(2*nm+nz)  ; 

F=sqrtm(L) ;Fi=F\eye(2*nm+nz) ; 

bl=bl*Fi ; cl=F*cl ; dll=F*dl l*Fi ; d21=d21*Fi ; dl2=F*dl2 ; 

y,  Form  the  outputs  needed  for  Gahinet's  klmi.m 
Ps=pck(a,[bl  b2] , [cl;c2] , [dll  dl2;d21  d22]); 

end 
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*/,  end  nshinflmi.m 


'/'/'/'/•/'/•/•/•/•/•/•/•/•/•/•/•/•/•/•/vvvvvvvyvvvvvvvv»/v»/v»/»/vvv»/vvv»/vv'/*/*/vv»i 


C.IO  nshinf  lmi2  .m;  Solve  LMI  GEVP  Problem 


iyxayMyx/x/x/x/x/x/x/x/x/x/x/xm 


'/,  nshinf lmi2.ni 

f unct ion  [Ps , R , S , L , J , gamma] =nshinf lmi2 (P , blk , nw , opt , target ) 

•/. 

%  NSHINFLMI2  Find  the  solution  to  the  Apkarian's  scaled  Hinf  optimization 


'/•  problem.  The  problem  is  formulated  as  a  Generalized  Eigenvalue 

%  Problem  using  LMIs.  Gamma  is  displayed  during  the 

'/•  minimization. 

•/. 

'/,  usage:  [R,S,L, J]=nshinflmi2(P, blk, nw, opt, target) 

•/. 

y,  inputs:  P  -  Open-loop  plant  in  packed  notation,  i.e.  P=pck(A,  ... 

'/.  [Bt  B1  B2],  [Ct;Cl;C2],  [Dtt  ...;Dlt  .  .  .  ;  D2t  ...]), 

'/•  where  t  denotes  the  parameter  influence 

y«  blk  -  nm  X  2  matrix  describing  the  structure  of  the  parameter 

'/•  block.  The  first  element  in  each  row  denotes  the  size  of 

%  the  subblock.  The  second  element  in  each  row  is  either 

'/•  0  (for  repeated  scalar)  or  1  (for  full)  . 

%  nw  -  Number  of  exogeneous  inputs  (must  equal  the  number  of 

'/•  controlled  outputs) 

*/•  opt  -  Optional  vector  of  parameter  to  pass  to  the  LMI 

'/•  generalized  eigenvalue  solver.  opt=[  (accuracy)  (iters) 

'/•  (feas.  rad.)  (slow  prog,  tol)  (print)],  def ault=  [le-2 

'/•  200  le7  30  0]  .  See  the  LMI  Lab  manual  for  more  details, 

y*  target-  Optional  target  for  the  generalized  eigenvalue 

'/•  minimization  problem  (gevp)  ,  i.e.  the  desired  geimma. 

y. 

y,  outputs:  Ps  -  The  COMPLETE  weighted  plant  scaled  by  L  and  J 

%  R,  S  -  Solutions  to  the  gevp  problem 

*/•  L,  J  -  The  complete  scaling  matrices  derived  from  L3  and  J3 

y»  gamma  -  The  optimal  gamma  found 
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'/,  Written  by:  Capt  Mark  Spillman 

•/.  WL/FIGC-3  Bldg  146,  WPAFB 

•/. 

*/,  Modified  by:  Capt  Martin  Breton,  AFIT 
if  nargin<l 

dispC'  usage:  [R,S,L,J,gopt]=nshinflmi2(P,blk,nw, opt, target) ; ; 

dispC'  O; 

return 

end 

if  nargin<5;target=0;end 
if  nargin<4;opt=[le-2  200  le7  30  0];end 

y,  The  equations  described  below  are  derived  from  "A  Convex  Characterization 
y,  of  Gain-Scheduled  Hinf  Controllers" 

y,  Get  nm  and  change  blk  to  denote  the  structure  of  L3  and  J3 

nm=sum(blk( : ,1))  ; 

blk(:,2)=-(blk(:,2)-l); 

y,  Set  nz=nw,  a  current  restriction  of  the  code 
nz=nw ; 

y,  Unpack  P  euid  spilt  apart  the  elements  in  B,  C  and  D 
[A,B,C,D] =unpck(P) ; 

na=size(A, 1) ;nu=size(B,2)-nw-nm;ny=size(C, l)-nz-nm; 

Bt=B( : , 1 :nm) ;B1=B( : ,nm+l :nm+nw) ;B2=B( : ,nm+nw+l :nm+nw+nu) ; 

Ct=C(l :nm, : ) ;Cl=C(nm+l :nm+nz, : ) ;C2=C(nm+nz+l :nm+nz+ny, : ) ; 

Dtt=D(l :nm, 1 :nm) ;Dtl=D(l :nm,nm+l :nm+nw) ;Dt2=D(l :nm,nm+nw+l :nm+nw+nu) ; 
Dlt=D(nm+l :nm+nz, 1 :nm) ;Dll=D(nm+l :nm+nz,nm+l :nm+nw) ; 

D12=D(nm+l :nm+nz,nm+nw+l :nm+nw+nu) ; 

D2t=D(nm+nz+l :nm+nz+ny, 1 :nm) ;D21=D(nm+nz+l :nm+nz+ny,nm+l :nm+nw) ; 
D22=D(nm+nz+l :nm+nz+ny,nm+nw+l :nm+nw+nu) ; 

y,  tolerances  (tedcen  from  Gahinet's  routines) 
macheps=mach_eps ; 
tolsing=sqrt (macheps) ; 
toleig=macheps“ (2/3) ; 

y,  Meike  the  appropriate  substitutions  to  use  a  part  of  Gahinet's 
y,  routine  below 
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DD12=[Dt2;D12] ;DD21=[D2t  D21] ; 


'/,  The  following  was  taken  from  Galiinet's  goptlmi.m: 

y.  For  numerical  stability  of  the  controller  computation, 

y,  zero  the  sing,  values  of  DD12  s.t  I  I  B2  DD12“+  |  |  >  1/tolsing 

[u,s,v]=svd(DD12) ; 

abstol=max(toleig*norm(B2,l) ,tolsing*s(l , 1)) ; 
ratio=max([s;zeros(l,size(s,2))]) ./. . . 

max( [tolsing*abs(B2*v) ;abstol*ones(l ,nu)] ) ; 
ind2=f ind(ratio  <  1);  12=length(ind2) ; 

if  12  >  0,  s( : ,ind2)=zeros(nz, length (ind2)) ;  DD12=u*s*v';  end 
[u , s , v] =svd (DD21 ' ) ; 

abstol=max(toleig*norm(C2,l) ,tolsing*s(l , 1)) ; 
ratio=max([s;zeros(l,size(s,2))]) ./. . . 

max([tolsing*abs(C2'*v) ;abstol*ones(l,ny)] ) ; 
ind2=f ind(ratio  <  1);  12=length(ind2) ; 

if  12  >  0,  s(: ,ind2)=zeros (nm+nw, length (ind2)) ;  DD21=v*s ' *u' ;  end 

^9|C3)Ct*3|C3tC>t<*>)t*4!’)C***!<:***i<<********3|:*********’l<>|C3(C****3|:**3|C3|‘***>l<********* 

y,  Compute  the  outer  factors  for  the  standard  feasibility  problem 
Nr=rnull([B2;DD12] ' ,0,tolsing) ; 

Wrl=Nr(l :na, :) ;Wr2=Nr(na+l :na+nm+nz, :) ; 

Ns=rnull([C2,DD21] ,0,tolsing) ; 

Wsl=Ns(l :na, :) ;Ws2=Ns(na+l:na+nm+nz, :) ; 

y,  Cut-down  Wr2  for  the  GEVP 
[Wr2u,Wr2s,Wr2v] =svd(Wr2) ; 

[rwr2,cwr2]=size(Wr2) ;Wr2sp=zeros(rwr2,cwr2) ; 

cut  J=sum(diag(Wr2s)  >max  (size  (Wr2)  )  *mcix  (diag  (Wr2s)  )  *eps)  ; 

Wr2  sp (rwr2-cut  J+ 1 : rwr2 , cwr2- cut  J+ 1 : cwr2 ) =Wr2  s ( 1 ; cut  J , 1 : cut  J ) ; 
Wr2up= [Wr2u( : ,cut J+1 :rwr2)  Wr2u( : , 1 : cut J)] ; 

Wr2vp=[Wr2v( : , cut J+1 :cwr2)  Wr2v( : , 1 : cut J)] ; 

Wr2p=Wr2up*Wr2sp ; 

Wr2pp=Wr2u(: ,l:cutJ)*Wr2s(l:cutJ,l:cutJ) ; 

y.  Cut-down  Ws2  for  the  GEVP 
[Ws2u,Ws2s,Ws2v] =svd(Ws2) ; 

[rws2,cws2]=size(Ws2) ;Ws2sp=zeros(rws2,cws2) ; 
cutL=sum(diag(Ws2s)>max (size (Ws2))*max (diag (Ws2s))*eps) ; 
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Ws2sp(rws2-cutL+l :rws2,cws2-cutL+l:cws2)=Ws2s(l :cutL,l :cutL) ; 
Ws2iip=[Ws2u(  :  ,cutL+l  :rws2)  Ws2u(:  ,l:cutL)]  ; 

Ws2vp=[Ws2v( : ,cutL+l :cws2)  Ws2v(: ,l:cutL)] ; 

Ws2p=Ws2up*Ws2sp ; 

Ws2pp=Ws2u(: ,l:cutL)*Ws2s(l:cutL,l:cutL) ; 

'/,  Compute  the  New  outer  factors  for  the  GEVP 
cnr=size(Nr,2) ; 

Nrp=[Wr2vp  zeros (cnr, nm+nz) ; zeros (nm+nz,cnr)  eye(nm+nz)] ; 
cns=size(Ns,2) ; 

Nsp=[Ws2vp  zeros (cns, nm+nz) ; zeros (nm+nz, cns)  eye (nm+nz) ] ; 

'/,  Define  some  additional  shortcut  matrices 
Blh=[Bt  Bl] ;Clh=CCt;Cl] ;Dllh=[Dtt  Dtl;Dlt  Dll]; 

Il=[eye(nm)  zeros(nm,nz)] ; 

12= [zeros (nm , nm+nz) ; zeros (nz ,nm)  eye (nz) ] ; 

IY1= [zeros (cut J,cwr2-cutJ)  eye (cut J)] ; 

IZ1= [zeros (cutL,cws2-cutL)  eye(cutL)] ; 

y,  Define  the  matrix  variables:  R,  S,  L3,  J3 
setlmis(  []) ; 

R=lmivar(l, [na  1]  ) ; 

S=lmivar(l, [na  1]); 

L3=lmivar(l ,blk) ; 

[J3,nstandardvar] =lmivar(l,blk) ; 

y,  Define  the  variable  Y 
[Yl,ndec]=lmivar(l, [cutJ  1]); 
dy2=ndec+l: (nm+nz) *cutJ+ndec; 
dy2= (reshape (dy2 , cut J , (nm+nz) ) ) ' ; 

Y2=lmivar(3, [zeros (nm+nz, cwr2-cutJ)  dy2] ) ; 

Y2p=lmivar(3,dy2) ; 

Y3=lmivar(l , [(nm+nz)  1]  ) ; 

y,  Define  the  variable  Z 
[Zl,ndec]=lmivar(l, [cutL  1]); 
dz2=ndec+l : (nm+nz) *cutL+ndec ; 
dz2= (reshape (dz2 , cutL , (nm+nz) ) ) ' ; 

Z2=lmivar(3, [zeros (nm+nz, cws2-cutL)  dz2] ) ; 

Z2p=lmivar(3,dz2) ; 

[Z3,ntotalvar]=lmivar(l, [(nm+nz)  1]) ; 
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'/,  Display  the  variable  information 
nduinmyvar=ntotalvar-nstandardvar ; 
dispC' 

dispC 'Attempting  to  solve  the  GEVP.') 

dispCE 'There  are  '  int2str(nstandeu:dvar)  '  staindard  variables .'] ) 
dispCE'There  are  '  int2str(ndummyvar)  '  additional  dummy  variables.']) 
dispC'  '); 

'/,  Define  individual  terms  in  the  first  equation 
lmiterm([-l,l,l,R]  ,1,1) ; 

ImitermC [-1 ,2, 1,0]  , 1) ; 

ImitermC [-1 ,2,2,S]  ,1,1); 

y,  Define  the  individual  terms  in  the  second  equation 
ImitermC [-2, 1,1, L3] ,1,1); 

ImitermC [-2, 2, 1,0] ,1) ; 

ImitermC [-2,2,2, J3] ,1,1); 

y,  Define  the  individual  terms  in  the  third  equation 
ImitermC [3, 0,0,0] ,Nrp) ; 
lmitermC[3,l,l,R] ,Hrl'*A,Wrl , 's' ) ; 

ImitermC [3, 1 , 1 ,R] ,Wr2 ' *Clh, Wrl , ' s ' ) ; 

Imit  erm  C[3,2,l,J3],Il',Il=t'CBlh'  *Wrl+Dl  Ih '  *Wr2)  )  ; 

Imit erm  C [3 , 2 , 1 , 0] , 12*  CBlh ' *Wrl+Dl Ih ' *Wr2) ) ; 
lmitermC[-3,l,l,Yl] ,IY1' ,IY1); 

ImitermC [-3,2, 1,Y2] ,1,1); 

ImitermC [-3,2,2,Y3] ,1,1); 

y,  Define  the  individual  terms  in  the  fourth  equation 
ImitermC [4, 0,0,0] ,Nsp) ; 

ImitermC [4,1 ,1,S] ,Wsl'*A' ,Wsl, 's') ; 

ImitermC [4, 1,1,S] ,Ws2'*Blh' ,Wsl, 's') ; 

ImitermC [4, 2, 1,L3] ,11' ,Il*CClh*Wsl+Dllh*Ws2)) ; 

ImitermC [4, 2, 1,0] ,I2*CClh*Wsl+Dllh*Ws2)) ; 

ImitermC [-4, 1,1, Zl] ,IZ1' ,IZ1); 

ImitermC [-4,2, 1,Z2] ,1,1); 

ImitermC [-4 ,2,2, Z3] ,1,1); 

y,  Define  the  individual  terms  in  the  fifth  equation 
ImitermC  [5, 1,1,Y1], 1,1); 

ImitermC [5,2 ,1 ,Y2p] ,1,1); 
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ImitermC [5,2 ,2, Y3] ,1,1); 
lmiteriii([-5,l,l,  J3]  ,Wr2pp'*Il'  ,Il*Wr2pp)  ; 

ImitermC  [-5, 1, 1,0] ,Wr2pp’ *I2*Wr2pp) ; 

ImitermC [-5,2,2, J3] ,11' ,11) ; 

ImitermC [-5, 2, 2,0] ,12) ; 

'/,  Define  the  individual  terms  in  the  sixth  equation 
ImitermC [6, 1 , 1,Z1] , 1 ,1)  ; 

ImitermC [6, 2, l,Z2p] ,1,1) ; 

ImitermC  [6,2 ,2, Z3] ,1,1); 

ImitermC [-6, 1, 1,L3] ,Ws2pp'*Il' ,Il*Ws2pp) ; 

ImitermC [-6, 1,1,0] ,Ws2pp' *I2*Ws2pp) ; 

ImitermC [-6,2,2,L3] , II ’ ,11) ; 

ImitermC  [-6, 2, 2,0] ,12) ; 

lgevp=getlmis; 

'/,  Solve  the  Generalized  Eigenvalue  Minimization  Problem 
[geutima,xopt]=gevpClgevp,2,opt,  []  ,  []  , target)  ; 

R=dec2mat  Clgevp , xopt , R) ; 

S=dec2mat  Clgevp , xopt , S) ; 

L3=dec2mat Clgevp, xopt, L3) ; 

J3=dec2mat Clgevp, xopt, J3) ; 

'/,  Compute  L  and  J  Cold  way) 

'/,N=  C  J3\ey  e  Cnm)  )  -L3 ; 

'/,L=[-N'  M  zeros  Cnm,  nz)  ;N'  L3  zerosCnm,nz)  ;  zeros  Cnz,2*nm)  eyeCnz)]  ; 
'/,J1=C"N’  )\eyeCnm)+J3; 

'/,J=[  J1  J3  zerosCnm,nz)  ;  J3  J3  zeros  Cnm,  nz)  ;  zeros  Cnz,2*nm)  eyeCnz)]  ; 

’/,  Compute  L  and  J  Cnew  way) 

J3p=cholCJ3) ; 

J3pi=J3p\eyeCnm) ; 

J3i=J3pi*J3pi' ; 

N=CJ3i-L3); 

Np=cholC-N) ; 

Npi=Np\eyeCnm) ; 

Ni=Npi*Npi' ; 

Jl=J3+Ni; 

L=[-N  N  zeros Cnm, nz) ;N  L3  zerosCnm,nz) ;  zerosCnz,2*nm)  eyeCnz)]; 

J=[  J1  J3  zerosCnm,nz) ;  J3  J3  zeros Cnm, nz) ; zeros Cnz,2*nm)  eyeCnz)]; 
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'/,  Form  the  elements  of  the  true  weighted  plant 
a=A;bl=[zeros(na,nm)  Bt  Bl];b2=[B2  zeros (na,nm)] ; 
cl=[zeros(nm,na) ;Ct;Cl] ;c2=[C2; zeros (nm,na)] ; 
dll= [zeros(nm,2*nm+nw) ;zeros(nm)  Dtt  Dtl ; zeros (nz,nm)  Dlt 
dl2= [zeros (nm,nu)  eye(nm);Dt2  zeros(nm) ;D12  zeros(nz,nm)] 
d21= [zeros(ny,nm)  D2t  D21;  eye(nm)  zeros (nm,nm+nw)] ; 
d22=[D22  zeros(ny,nm) ; zeros (nm,nu+nm)] ; 

'/,  Scale  the  true  plant  so  that  L=J=I 

condL=cond(L) 

condJ=cond(J) 

F=sqrtm(L) ;Fi=F\eye(2*nm+nz) ; 

'/,F=chol(L)  ;Fi=F\eye(2*nm+nz) ; 

bl=bl*Fi ; cl=F*cl ; dl l=F*dll*Fi ; d21=d21*Fi ; dl2=F*dl2 ; 

'/,  Form  the  outputs  needed  for  Gahinet's  klmi.m 
Ps=pck(a,[bl  b2] , [cl;c2] , [dll  dl2;d21  d22]); 


'/•  end  nshinflmi2.m 


Appendix  D.  Simulation  Programs 


This  appendix  lists  the  main  programs  and  S-functions  used  to  setup  and  run 
the  simulations. 


D.l  LPVG.m;  Generate  LPV  Simulation  Aircraft  Model  in  SIMULINK 


•/.  LPVG.m 

function  [sys,xO]=LPVG(t ,x,u, flag, G,Mbarl,Mt ill, hbarl,htill,nMl ,nhl,Mo,ho,Uo) 

'/,  M-file  Simulink  S-function  to  incorporate 
y,  theta  into  G  before  converting  to  state-space 

'/,  nominal  G  must  be  in  workspace 

if  flag  ==  0 
M=Mo; 
h=ho; 
else 
M=u(l); 

h=u(2) ;  y  new  ho;  could  find  new  Uo 
uu=[u(3) ;u(4)] ; 
end 

dM=(M-Mbarl)/Mtill; 

dh=(h-hbarl)/htill; 

theta=[dM*eye(nMl)  zeros (nMl ,nhl) ; zeros (nhl,nMl)  dh*eye(nhl)] ; 

GLTI=starp (theta, G) ; 

[Ag,Bg,Cg,Dg]=unpck(GLTI) ; 

y,  add  wind  noise  input  for  simulink; 

Bgx=[Ag(:,2)  Bg] ; 

Dgx= [zeros (4,1)  Dg] ; 

if  abs(flag)  ==  1  51  return  state  derivatives 
sys=Ag*x+Bgx*uu ; 
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elseif  flag  ==  3  */,  return  system  outputs  y 

sys=Cg*x+Dgx*uu;  %  y  must  be  at  least:  u,aoa, theta 
x(l)=0; 

elseif  flag  ==  0  '/,  parameter  sizes  and  initial  conditions 

[m,n]=size(Dgx) ; 

sys=  [length (Ag)  0  m  n+2  0  meix(any(Dgx~=0))]  ; 
x0=[0  0  0  0]^ 

else  y,  flag=2  or  4  used  for  discrete  only 
sys=[]  ; 
end 


•/.  end  of  LPVG.m 


D.2  LPVK.m;  Generate  LPV  Controller  in  SIMULINK 


y,  LPVK.m 

function  [sys ,xO]=LPVK(t ,x,u,flag,K,Mbar,Mtil,hbar ,htil ,nM,nh,Mo,ho) 

y,  M-file  Simulink  S-function  to  incorporate 
y,  theta  into  K  before  converting  to  state-space 
y,  use  in  dynamic  simulations 

y,  nominal  K  must  be  in  workspace 

if  flag  ==  0 
M=Mo; 
h=ho; 
else 
M=u(l) ; 
h=u(2) ; 

uu=[u(3) ;u(4) ;u(5)] ; 
end 
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dM=(M-Mbar)/Mtil; 
dh= (h-hbar) /htil ; 

theta=[dM*eye(nM)  zeros(nM,nh)  ;2eros(nh,iiM)  dh*eye(nh)]  ; 

KLTI=starp(K, theta) ; 

[Ak,Bk,Ck,Dk] =unpck(KLTI) ; 

if  abs(flag)  ==  1  '/,  return  state  derivatives 
sys=Ak*x+Bk*uu ; 

elseif  flag  ==  3  JC  return  system  outputs  y 

sys=Ck*x+Dk*uu ;  '/,  y  must  be  at  least:  u,aoa, theta 

elseif  flag  ==  0  '/,  parameter  sizes  and  initial  conditions 

[m,n]=size(Dk) ; 

sys=[length(Ak)  0  m  n+2  0  max(any(Dk“=0) )] ; 
xO=zeros (length (Ak) , 1) ; 

else  '/,  flag=2  or  4  used  for  discrete  only 
sys=[]  ; 
end 

•/.  end  of  LPVK.m 


nvixav‘ 


D.3  f  ISsimdat  .m;  Setup  Static  Simulations 


•/•/•/•/•/•/•/»/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/»/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/»/•/•/•> 


'/,  f ISsimdat. m 

y,  fl8  static  simulation  setup  program 

y,  this  program  converts  the  LPV  plant 
y,  eind  LPV  controller  into  an  LTI  plant 
y,  and  LTI  controller  at  a  given  M,h  pt 
y,  in  the  flight  envelope  in  order  to 
y,  perform  static  tests 
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y,  load  LPV  K  y,  need  controller  in  workspace  or  file 
load  fenvSRP 

y,  enter  given  values  for  M,h  test  point 
M=.75 
h=15000 
ao=0 

y,  end  of  user  inputs 

K=K_3; 
nM=tblk(l); 
nh=tblk(2) ; 

Act=[  -20.2000  20.2000 

1.0000  0 

20.2000  20.2000 

0  0 

[Ak , Bk , Ck , Dk] =unpck (K) ; 

[Ade , Bde , Cde , Dde] =unpck (Act ) ; 

y,  LTI  Gdel  cind  Kdel  (which  taike  theta  into  account) 

dM=(M-Mbar)/Mtil; 

dh= (h-hbar) /ht il ; 

theta=[dM*eye(nM)  zeros (nM, nh) ; zeros (nh,nM)  dh*eye(nh)] ; 

Kdel=starp(K, theta) ;  '/,  gives  K  at  test  point 
[Ak,Bk,Ck,Dk] =unpck(Kdel) ; 

y,  dele  actuator  output 
Cde=Cde(l,l) ; 

Dde=Dde(l,l); 


1.0000 

0 

0 

-Inf]  ; 


y,  get  full  longitudinal  simulation  model  and  find  LTI  G  at  test  pt 

load  LPVF18;  '/,  get  full  longitudinal  plant  G 

nMl=9; 

nhl=10; 

y,  LTI  Gdel  (which  take  theta  into  account) 
dMl=(M-Mbarl)/Mtill; 


D-4 


dhl=(h-hbarl)/htill ; 

thetl=[dMl*eye(nMl)  zeros (nMl, nhl) ; zeros (nhl ,nMl)  dhl*eye(nhl)] ; 

Gdel=starp(thetl ,G) ; 

[Ag ,  Bg , Cg , Dg]  =mipck ( Gdel)  ; 

s=size(Ag,l)/2;  '/,  s  column  of  Ag  is  aoa 
'/,  add  wind  noise  input  for  simulink; 

Bgx=[Ag(:,s)  Bg] ; 

Dgx=[zeros(size(Dg,l) ,1)  Dg] ; 

y,  end  of  flBsimdat.m 


y//////////.y;///////////.y//////////.y//////////.y// 


';/.y//////////////////////.y//////.y.y//////.y.y.y.y.y. 


D.4  flSLPVsimdat  .m;  Setup  Dynamic  Simulation 


y,  f  1 8LPV s  imdat .  m 

y,  fl8  dynamic  simulation  setup  program 

y,  this  program  sets  up  LPV  plant  simulation  model 

y,  in  order  to  perform  dynamic  tests 

y,  controller  must  be  available  in  workspace  or  file 

load  envleB 

K=K_3;  y,  LPV  controller 
'/.load  Kmusp 
'/.load  Kmuhih 

y,[Ak,Bk,Ck,Dk]=unpck(k_dk3gb) ;  '/.  to  test  LTI/mu  controller 

y,  need  to  also  replace  LPVK  with  state-space 

nM=tblk(l); 

nh=tblk(2) ; 

hbar 

htil 

Mbar 

Mtil 
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y,  user  inputs  —  starting  pt  of  flight  trajectory 

ho=30000 

Mo=.5 

ao=0 

y,  end  user  inputs 

To=518.67-.003565*ho  good  approx  up  to  36000  ft 
Uo=Mo*sqrt (1. 4*1716. 16*To) 

load  LPVF18;  '/,  get  full  longitudinal  simulation  plant 

nMl=9; 

nhl=10; 

Act=[  -20.2000  20.2000  1.0000 

1.0000  0  0 

20.2000  20.2000  0 

0  0  -Inf] ; 

[Ade , Bde , Cde , Dde] =unpck (Act ) ; 
y,  dele  actuator  output 
Cde=Cde(l,l) ; 

Dde=Dde(l,l) ; 

y,  end  of  f  18LPVsimdat  .m 
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Systems  which  vary  significantly  over  an  operating  envelope,  such  as  fighter  aircraft,  generally  cannot  be 
controlled  by  a  single  linear  time-invariant  controller.  As  a  result,  gain-scheduling  methods  are  employed 
to  design  control  laws  which  can  provide  the  desired  performance.  This  thesis  examines  a  relatively  new 
approach  to  gain-scheduling,  in  which  the  varying  controller  is  designed  from  the  outset  to  guarantee  robust 
performance,  thereby  avoiding  the  disadvantages  of  point  designs.  Specifically,  the  parameter-varying  (LPV) 
aircraft  model  is  linearized  using  linear  fractional  transformations  (LFT’s),  and  the  resulting  control  problem 
is  characterized  as  the  solution  to  a  set  of  four  linear  matrix  inequalities  (LMI’s).  The  supporting  theory  is 
reviewed  and  two  pitch-rate  controllers  are  designed;  one  for  the  full  longitudinal  aircraft  model,  and  another 
for  the  short  period  model.  It  is  found  that,  even  though  the  varying  controllers  are  quite  conservative,  they 
can  guarantee  better  robust  performance  over  a  large  portion  of  an  operating  envelope  when  compared  to 
time-invariant  //-synthesis  controllers. 
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